{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to atomman: LAMMPS functionality\n",
    "\n",
    "__Lucas M. Hale__, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), _Materials Science and Engineering Division, NIST_.\n",
    "    \n",
    "[Disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction<a id='section1'></a>\n",
    "\n",
    "This Notebook provides an introduction to interacting with LAMMPS using atomman by working through a demonstration simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Library Imports**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atomman version = 1.4.0\n",
      "Notebook executed on 2021-08-04\n"
     ]
    }
   ],
   "source": [
    "# Standard libraries\n",
    "from pathlib import Path\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np             \n",
    "\n",
    "# http://matplotlib.org/\n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline\n",
    "\n",
    "# https://github.com/usnistgov/atomman\n",
    "import atomman as am            \n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# Show atomman version\n",
    "print('atomman version =', am.__version__)\n",
    "\n",
    "# Show date of Notebook execution\n",
    "print('Notebook executed on', datetime.date.today())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Interatomic Potential Control<a id='section2'></a>\n",
    "\n",
    "LAMMPS is capable of using a wide array of interatomic potential styles, which are defined through a combination of LAMMPS commands and potential parameter files.  In atomman, the LAMMPS commands can be automatically generated using the Potentials class and structured data model files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Load potential from database\n",
    "\n",
    "Any LAMMPS-compatible potential in the NIST Interatomic Potentials Repository can be loaded/downloaded using load_lammps_potential().  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load potential based on its unique ID\n",
    "pot_id = '2015--Pascuet-M-I--Al--LAMMPS--ipr1'\n",
    "potential = am.load_lammps_potential(id=pot_id, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Potential parameter files\n",
    "\n",
    "Setting getfiles=True will copy/download the potential parameter files from the local/remote library database to the current working directory in a subfolder matching the potential's id."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2015--Pascuet-M-I--Al--LAMMPS--ipr1\\Al.meam\n",
      "re(1,1) =   2.86378\n",
      "attrac(1,1) = 0.39558\n",
      "repuls(1,1) = 0.09806\n",
      "Cmin(1,1,1) = 1.00769\n",
      "Cmax(1,1,1) = 2.31407\n",
      "#\n",
      "rc = 5.0\n",
      "delr = 0.1\n",
      "augt1 = 1\n",
      "\n",
      "\n",
      "2015--Pascuet-M-I--Al--LAMMPS--ipr1\\library-Al.meam\n",
      "# References:\n",
      "# elt        lat     z       ielement     atwt\n",
      "# alpha      b0      b1      b2           b3    alat    esub    asub\n",
      "# t0         t1              t2           t3            rozero  ibar\n",
      "#\n",
      "'Al'        'fcc'   12.     13           26.9815\n",
      "  4.68604   1.56205   5.39270   5.29601  -1.00047   4.05000   3.35999   1.06859\n",
      "1.0  -1.54917  -1.28508  10.01041 1.0 0 \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for fname in Path(pot_id).glob('*'):\n",
    "    print(fname)\n",
    "    with open(fname) as f:\n",
    "        print(f.read())\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. lammps.Potential Class\n",
    "\n",
    "The load_lammps_potential() function returns a Potential object that can be used to explore properties of the potential and generate LAMMPS input commands."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic properties can be directly obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "str(potential) ->        potential_LAMMPS record named 2015--Pascuet-M-I--Al--LAMMPS--ipr1\n",
      "potential.units ->       metal\n",
      "potential.atom_style ->  atomic\n",
      "potential.pair_style ->  meam\n",
      "potential.symbols ->     ['Al']\n",
      "potential.elements() ->  ['Al']\n",
      "potential.masses() ->    [26.9815]\n"
     ]
    }
   ],
   "source": [
    "# Show that basic properties can be directly obtained\n",
    "print('str(potential) ->       ', potential)\n",
    "print('potential.units ->      ', potential.units)\n",
    "print('potential.atom_style -> ', potential.atom_style)\n",
    "print('potential.pair_style -> ', potential.pair_style)\n",
    "print('potential.symbols ->    ', potential.symbols)\n",
    "print('potential.elements() -> ', potential.elements())\n",
    "print('potential.masses() ->   ', potential.masses())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LAMMPS command lines for the potential are auto-generated based on a list of symbols corresponding to LAMMPS atom types. This fully works for all LAMMPS pair_styles, with only the hybrid and hybrid/overlay styles having limitations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "potential.pair_info() ->\n",
      "print \"Potential 2015--Pascuet-M-I--Al--LAMMPS--ipr1 listed in the NIST Interatomic Potentials Repository:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/entry/2015--Pascuet-M-I-Fernandez-J-R--Al/2015--Pascuet-M-I--Al--LAMMPS--ipr1.html\"\n",
      "print \"Publication(s) related to the potential:\"\n",
      "print \"https://doi.org/10.1016/j.jnucmat.2015.09.030\"\n",
      "print \"Parameter file(s) can be downloaded at:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/Al.meam\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/library-Al.meam\"\n",
      "pair_style meam\n",
      "pair_coeff * * 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\library-Al.meam Al 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\Al.meam Al\n",
      "mass 1 26.9815\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('potential.pair_info() ->')\n",
    "print(potential.pair_info())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "potential.pair_info(['Al', 'Al', 'Al']) ->\n",
      "print \"Potential 2015--Pascuet-M-I--Al--LAMMPS--ipr1 listed in the NIST Interatomic Potentials Repository:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/entry/2015--Pascuet-M-I-Fernandez-J-R--Al/2015--Pascuet-M-I--Al--LAMMPS--ipr1.html\"\n",
      "print \"Publication(s) related to the potential:\"\n",
      "print \"https://doi.org/10.1016/j.jnucmat.2015.09.030\"\n",
      "print \"Parameter file(s) can be downloaded at:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/Al.meam\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/library-Al.meam\"\n",
      "pair_style meam\n",
      "pair_coeff * * 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\library-Al.meam Al 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\Al.meam Al Al Al\n",
      "mass 1 26.9815\n",
      "mass 2 26.9815\n",
      "mass 3 26.9815\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\"potential.pair_info(['Al', 'Al', 'Al']) ->\")\n",
    "print(potential.pair_info(['Al', 'Al', 'Al']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Generate initial system<a id='section3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Load relaxed crystal \n",
    "\n",
    "A crystalline system can be easily generated using a unit cell system either defined in atomman, imported from another format, or accessed from the potentials database.  For simplicity, here we will get the relaxed fcc structure for the potential."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load fcc prototype with Al lattice constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 4.050,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  4.050,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  4.050]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Al',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   2.025 |   2.025\n",
      "      2 |       1 |   2.025 |   0.000 |   2.025\n",
      "      3 |       1 |   2.025 |   2.025 |   0.000\n"
     ]
    }
   ],
   "source": [
    "system = am.load('crystal', family='A1--Cu--fcc', potential=potential)\n",
    "print(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Manipulate system\n",
    "\n",
    "More complicated atomic configurations can then be generated by manipulating the seed system and the atoms contained within. Here, we'll limit the manipulations to making the system a 3x3x3 supercell of itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Supercell has 108 atoms\n"
     ]
    }
   ],
   "source": [
    "system = system.supersize(3,3,3)\n",
    "print('Supercell has', system.natoms, 'atoms')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Save to atom data file\n",
    "\n",
    "System.dump('atom_data') outputs the system to a LAMMPS atom data file. Quick notes on the parameters used here\n",
    "\n",
    "- Giving potential allows for the appropriate units and atom_style settings to be used.\n",
    "- The float_format value used here is a small precision to enhance clarity of print statements below. Typically, the precision should be large (default value is '%.13f').\n",
    "- Setting return_pair_info=True will return the LAMMPS commands for both the system and the potential together.  This is the preferred method as it ensures compatibility with all potential pair styles currently in the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save System to 'atom.dat' atom data file \n",
    "system_pair_info = system.dump('atom_data', f='atom.dat',\n",
    "                          potential=potential,\n",
    "                          float_format='%.4f',    # Remove or make larger precision for real runs!\n",
    "                          )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the returned LAMMPS command lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "units metal\n",
      "atom_style atomic\n",
      "\n",
      "boundary p p p\n",
      "read_data atom.dat\n",
      "\n",
      "print \"Potential 2015--Pascuet-M-I--Al--LAMMPS--ipr1 listed in the NIST Interatomic Potentials Repository:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/entry/2015--Pascuet-M-I-Fernandez-J-R--Al/2015--Pascuet-M-I--Al--LAMMPS--ipr1.html\"\n",
      "print \"Publication(s) related to the potential:\"\n",
      "print \"https://doi.org/10.1016/j.jnucmat.2015.09.030\"\n",
      "print \"Parameter file(s) can be downloaded at:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/Al.meam\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/library-Al.meam\"\n",
      "pair_style meam\n",
      "pair_coeff * * 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\library-Al.meam Al 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\Al.meam Al\n",
      "mass 1 26.9815\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(system_pair_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the contents of the data file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "108 atoms\n",
      "1 atom types\n",
      "0.0000 12.1500 xlo xhi\n",
      "0.0000 12.1500 ylo yhi\n",
      "0.0000 12.1500 zlo zhi\n",
      "\n",
      "Atoms # atomic\n",
      "\n",
      "1 1 0.0000 0.0000 0.0000\n",
      "2 1 0.0000 2.0250 2.0250\n",
      "3 1 2.0250 0.0000 2.0250\n",
      "4 1 2.0250 2.0250 0.0000\n",
      "5 1 4.0500 0.0000 0.0000\n",
      "6 1 4.0500 2.0250 2.0250\n",
      "7 1 6.0750 0.0000 2.0250\n",
      "8 1 6.0750 2.0250 0.0000\n",
      "9 1 8.1000 0.0000 0.0000\n",
      "10 1 8.1000 2.0250 2.0250\n",
      "11 1 10.1250 0.0000 2.0250\n",
      "12 1 10.1250 2.0250 0.0000\n",
      "13 1 0.0000 4.0500 0.0000\n",
      "14 1 0.0000 6.0750 2.0250\n",
      "15 1 2.0250 4.0500 2.0250\n",
      "16 1 2.0250 6.0750 0.0000\n",
      "17 1 4.0500 4.0500 0.0000\n",
      "18 1 4.0500 6.0750 2.0250\n",
      "19 1 6.0750 4.0500 2.0250\n",
      "20 1 6.0750 6.0750 0.0000\n",
      "21 1 8.1000 4.0500 0.0000\n",
      "22 1 8.1000 6.0750 2.0250\n",
      "23 1 10.1250 4.0500 2.0250\n",
      "24 1 10.1250 6.0750 0.0000\n",
      "25 1 0.0000 8.1000 0.0000\n",
      "26 1 0.0000 10.1250 2.0250\n",
      "27 1 2.0250 8.1000 2.0250\n",
      "28 1 2.0250 10.1250 0.0000\n",
      "29 1 4.0500 8.1000 0.0000\n",
      "30 1 4.0500 10.1250 2.0250\n",
      "31 1 6.0750 8.1000 2.0250\n",
      "32 1 6.0750 10.1250 0.0000\n",
      "33 1 8.1000 8.1000 0.0000\n",
      "34 1 8.1000 10.1250 2.0250\n",
      "35 1 10.1250 8.1000 2.0250\n",
      "36 1 10.1250 10.1250 0.0000\n",
      "37 1 0.0000 0.0000 4.0500\n",
      "38 1 0.0000 2.0250 6.0750\n",
      "39 1 2.0250 0.0000 6.0750\n",
      "40 1 2.0250 2.0250 4.0500\n",
      "41 1 4.0500 0.0000 4.0500\n",
      "42 1 4.0500 2.0250 6.0750\n",
      "43 1 6.0750 0.0000 6.0750\n",
      "44 1 6.0750 2.0250 4.0500\n",
      "45 1 8.1000 0.0000 4.0500\n",
      "46 1 8.1000 2.0250 6.0750\n",
      "47 1 10.1250 0.0000 6.0750\n",
      "48 1 10.1250 2.0250 4.0500\n",
      "49 1 0.0000 4.0500 4.0500\n",
      "50 1 0.0000 6.0750 6.0750\n",
      "51 1 2.0250 4.0500 6.0750\n",
      "52 1 2.0250 6.0750 4.0500\n",
      "53 1 4.0500 4.0500 4.0500\n",
      "54 1 4.0500 6.0750 6.0750\n",
      "55 1 6.0750 4.0500 6.0750\n",
      "56 1 6.0750 6.0750 4.0500\n",
      "57 1 8.1000 4.0500 4.0500\n",
      "58 1 8.1000 6.0750 6.0750\n",
      "59 1 10.1250 4.0500 6.0750\n",
      "60 1 10.1250 6.0750 4.0500\n",
      "61 1 0.0000 8.1000 4.0500\n",
      "62 1 0.0000 10.1250 6.0750\n",
      "63 1 2.0250 8.1000 6.0750\n",
      "64 1 2.0250 10.1250 4.0500\n",
      "65 1 4.0500 8.1000 4.0500\n",
      "66 1 4.0500 10.1250 6.0750\n",
      "67 1 6.0750 8.1000 6.0750\n",
      "68 1 6.0750 10.1250 4.0500\n",
      "69 1 8.1000 8.1000 4.0500\n",
      "70 1 8.1000 10.1250 6.0750\n",
      "71 1 10.1250 8.1000 6.0750\n",
      "72 1 10.1250 10.1250 4.0500\n",
      "73 1 0.0000 0.0000 8.1000\n",
      "74 1 0.0000 2.0250 10.1250\n",
      "75 1 2.0250 0.0000 10.1250\n",
      "76 1 2.0250 2.0250 8.1000\n",
      "77 1 4.0500 0.0000 8.1000\n",
      "78 1 4.0500 2.0250 10.1250\n",
      "79 1 6.0750 0.0000 10.1250\n",
      "80 1 6.0750 2.0250 8.1000\n",
      "81 1 8.1000 0.0000 8.1000\n",
      "82 1 8.1000 2.0250 10.1250\n",
      "83 1 10.1250 0.0000 10.1250\n",
      "84 1 10.1250 2.0250 8.1000\n",
      "85 1 0.0000 4.0500 8.1000\n",
      "86 1 0.0000 6.0750 10.1250\n",
      "87 1 2.0250 4.0500 10.1250\n",
      "88 1 2.0250 6.0750 8.1000\n",
      "89 1 4.0500 4.0500 8.1000\n",
      "90 1 4.0500 6.0750 10.1250\n",
      "91 1 6.0750 4.0500 10.1250\n",
      "92 1 6.0750 6.0750 8.1000\n",
      "93 1 8.1000 4.0500 8.1000\n",
      "94 1 8.1000 6.0750 10.1250\n",
      "95 1 10.1250 4.0500 10.1250\n",
      "96 1 10.1250 6.0750 8.1000\n",
      "97 1 0.0000 8.1000 8.1000\n",
      "98 1 0.0000 10.1250 10.1250\n",
      "99 1 2.0250 8.1000 10.1250\n",
      "100 1 2.0250 10.1250 8.1000\n",
      "101 1 4.0500 8.1000 8.1000\n",
      "102 1 4.0500 10.1250 10.1250\n",
      "103 1 6.0750 8.1000 10.1250\n",
      "104 1 6.0750 10.1250 8.1000\n",
      "105 1 8.1000 8.1000 8.1000\n",
      "106 1 8.1000 10.1250 10.1250\n",
      "107 1 10.1250 8.1000 10.1250\n",
      "108 1 10.1250 10.1250 8.1000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "with open('atom.dat') as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Converting to/from LAMMPS units<a id='section4'></a>\n",
    "\n",
    "LAMMPS performs its calculations with values in one of multiple sets of pre-defined units.  The atomman.lammps.style submodule has some useful functions when working with different units options.\n",
    "\n",
    "**atomman.lammps.style.timestep()**\n",
    "\n",
    "The lammps.style.timestep() function returns the default timestep value for a given LAMMPS units option.\n",
    "\n",
    "Parameters\n",
    "\n",
    "- **units** (*str*) the LAMMPS units option being used.\n",
    "\n",
    "**atomman.lammps.style.unit()**\n",
    "\n",
    "The lammps.style.unit() function returns a dictionary giving the units associated with physical quantities as used by LAMMPS with a given units option.\n",
    "\n",
    "Parameters\n",
    "\n",
    "- **units** (*str*) the LAMMPS units option being used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.001\n"
     ]
    }
   ],
   "source": [
    "timestep = lmp.style.timestep(potential.units)\n",
    "print(timestep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OrderedDict([('mass', 'g/mol'), ('length', 'angstrom'), ('time', 'ps'), ('energy', 'eV'), ('velocity', 'angstrom/ps'), ('force', 'eV/angstrom'), ('torque', 'eV'), ('temperature', 'K'), ('pressure', 'bar'), ('dynamic viscosity', 'Pa*s/10'), ('charge', 'e'), ('dipole', 'e*angstrom'), ('electric field', 'V/angstrom'), ('density', 'g/cm^3'), ('ang-mom', 'angstrom*angstrom/ps*g/mol'), ('ang-vel', '1/ps')])\n"
     ]
    }
   ],
   "source": [
    "lammps_unit = lmp.style.unit(potential.units)\n",
    "print(lammps_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Composing LAMMPS Input Scripts<a id='section5'></a>\n",
    "\n",
    "LAMMPS scripts can easily be constructed by combining the system_info generated from System.dump('atom_data'), the pair_info from Potential.pair_info(), and any user-defined input lines.  This allows for specific simulation actions to easily be perfored across different potentials or initial configurations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. Show content generated in previous sections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show system_pair_info generated in Section 3.3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "units metal\n",
      "atom_style atomic\n",
      "\n",
      "boundary p p p\n",
      "read_data atom.dat\n",
      "\n",
      "print \"Potential 2015--Pascuet-M-I--Al--LAMMPS--ipr1 listed in the NIST Interatomic Potentials Repository:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/entry/2015--Pascuet-M-I-Fernandez-J-R--Al/2015--Pascuet-M-I--Al--LAMMPS--ipr1.html\"\n",
      "print \"Publication(s) related to the potential:\"\n",
      "print \"https://doi.org/10.1016/j.jnucmat.2015.09.030\"\n",
      "print \"Parameter file(s) can be downloaded at:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/Al.meam\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/library-Al.meam\"\n",
      "pair_style meam\n",
      "pair_coeff * * 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\library-Al.meam Al 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\Al.meam Al\n",
      "mass 1 26.9815\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(system_pair_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2. Write LAMMPS input script template\n",
    "\n",
    "LAMMPS scripts can be dynamically generated using Python functions or templates that take the above info lines and other values as parameters. Here, we demonstrate a LAMMPS input template script where all fields to be filled in in Python are delimited with <> brackets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "template = \"\"\"\n",
    "<system_pair_info>\n",
    "\n",
    "# Define temperature and dependent variables\n",
    "variable T equal <temperature>\n",
    "variable twoT equal 2*$T\n",
    "\n",
    "# Define equilibrium pressure\n",
    "variable P equal <pressure>\n",
    "\n",
    "# Define timestep and dependent variables\n",
    "variable deltat equal <timestep>\n",
    "variable Trelax equal 100*${deltat}\n",
    "variable Prelax equal 1000*${deltat}\n",
    "\n",
    "# Initialize atomic velocities with twoT\n",
    "velocity all create ${twoT} 124352\n",
    "\n",
    "# Define thermo steps and properties\n",
    "thermo 100\n",
    "thermo_style custom step temp press lx ly lz\n",
    "\n",
    "# Define dump \n",
    "dump mydump all atom 100000 *.dump\n",
    "\n",
    "# Specify timestep\n",
    "timestep ${deltat}\n",
    "\n",
    "# Apply npt conditions\n",
    "fix 1 all npt temp $T $T ${Trelax} iso $P $P ${Prelax}\n",
    "\n",
    "# Run simulation\n",
    "run 100000\n",
    "\"\"\"           "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3. Fill in the template with atomman.tools.filltemplate()\n",
    "\n",
    "The template can then be easily filled in with the atomman.tools.filltemplate() function. \n",
    "\n",
    "Parameters\n",
    "\n",
    "- **template** (*str or file-like object*) is the template to fill in.\n",
    "\n",
    "- **variable** (*dict*) gives the delimited keys and corresponding values to insert into the template.\n",
    "\n",
    "- **s_delimiter** (*str*) the starting delimiter for identifying variable names.\n",
    "\n",
    "- **e_delimiter** (*str*) the ending delimiter for identifying variable names."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build dictionary of template variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'system_pair_info': 'units metal\\natom_style atomic\\n\\nboundary p p p\\nread_data atom.dat\\n\\nprint \"Potential 2015--Pascuet-M-I--Al--LAMMPS--ipr1 listed in the NIST Interatomic Potentials Repository:\"\\nprint \"https://www.ctcms.nist.gov/potentials/entry/2015--Pascuet-M-I-Fernandez-J-R--Al/2015--Pascuet-M-I--Al--LAMMPS--ipr1.html\"\\nprint \"Publication(s) related to the potential:\"\\nprint \"https://doi.org/10.1016/j.jnucmat.2015.09.030\"\\nprint \"Parameter file(s) can be downloaded at:\"\\nprint \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/Al.meam\"\\nprint \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/library-Al.meam\"\\npair_style meam\\npair_coeff * * 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\\\library-Al.meam Al 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\\\Al.meam Al\\nmass 1 26.9815\\n\\n', 'timestep': 0.001, 'temperature': 100, 'pressure': 0.0}\n"
     ]
    }
   ],
   "source": [
    "lammps_variable = {}\n",
    "\n",
    "# Generated above\n",
    "lammps_variable['system_pair_info'] = system_pair_info\n",
    "\n",
    "# Set timestep to default value for LAMMPS units\n",
    "lammps_variable['timestep'] = lmp.style.timestep(units=potential.units)\n",
    "\n",
    "# Specify temperature to equilibriate at (always in Kelvin)\n",
    "lammps_variable['temperature'] = 100\n",
    "\n",
    "# Specify pressure to equilibriate at\n",
    "pressure = uc.set_in_units(0.0, 'MPa')\n",
    "lammps_variable['pressure'] = uc.get_in_units(pressure, lammps_unit['pressure'])\n",
    "\n",
    "print(lammps_variable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill in template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "units metal\n",
      "atom_style atomic\n",
      "\n",
      "boundary p p p\n",
      "read_data atom.dat\n",
      "\n",
      "print \"Potential 2015--Pascuet-M-I--Al--LAMMPS--ipr1 listed in the NIST Interatomic Potentials Repository:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/entry/2015--Pascuet-M-I-Fernandez-J-R--Al/2015--Pascuet-M-I--Al--LAMMPS--ipr1.html\"\n",
      "print \"Publication(s) related to the potential:\"\n",
      "print \"https://doi.org/10.1016/j.jnucmat.2015.09.030\"\n",
      "print \"Parameter file(s) can be downloaded at:\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/Al.meam\"\n",
      "print \"https://www.ctcms.nist.gov/potentials/Download/2015--Pascuet-M-I-Fernandez-J-R--Al/1/library-Al.meam\"\n",
      "pair_style meam\n",
      "pair_coeff * * 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\library-Al.meam Al 2015--Pascuet-M-I--Al--LAMMPS--ipr1\\Al.meam Al\n",
      "mass 1 26.9815\n",
      "\n",
      "\n",
      "\n",
      "# Define temperature and dependent variables\n",
      "variable T equal 100\n",
      "variable twoT equal 2*$T\n",
      "\n",
      "# Define equilibrium pressure\n",
      "variable P equal 0.0\n",
      "\n",
      "# Define timestep and dependent variables\n",
      "variable deltat equal 0.001\n",
      "variable Trelax equal 100*${deltat}\n",
      "variable Prelax equal 1000*${deltat}\n",
      "\n",
      "# Initialize atomic velocities with twoT\n",
      "velocity all create ${twoT} 124352\n",
      "\n",
      "# Define thermo steps and properties\n",
      "thermo 100\n",
      "thermo_style custom step temp press lx ly lz\n",
      "\n",
      "# Define dump \n",
      "dump mydump all atom 100000 *.dump\n",
      "\n",
      "# Specify timestep\n",
      "timestep ${deltat}\n",
      "\n",
      "# Apply npt conditions\n",
      "fix 1 all npt temp $T $T ${Trelax} iso $P $P ${Prelax}\n",
      "\n",
      "# Run simulation\n",
      "run 100000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Generate script from template and lammps_variable\n",
    "script = am.tools.filltemplate(template, lammps_variable, '<', '>')\n",
    "\n",
    "# Save script to 'nvt.in'\n",
    "with open('nvt.in', 'w') as f:\n",
    "    f.write(script)\n",
    "\n",
    "# Show contents of 'nvt.in'    \n",
    "with open('nvt.in') as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Run LAMMPS<a id='section6'></a>\n",
    "\n",
    "The LAMMPS simulation can be ran from within Python using the run() function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify your lammps executable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_exe = 'lmp_serial'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simulation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = lmp.run(lammps_exe, 'nvt.in')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting simulation data is returned as a Log object, which containes the thermo data from the log.lammps files. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Analyzing Thermo Results<a id='section7'></a>\n",
    "\n",
    "Data for each simulation run/minimization is stored in the returned Log object.  Each Simulation has a thermo property that is a pandas.DataFrame of the LAMMPS thermo data.\n",
    "\n",
    "*Updated version 1.3.7*: Each simulation is now represented using a Simulation class rather than a dictionary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show thermo data column names for the first (and only) simulation run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Step', 'Temp', 'Press', 'Lx', 'Ly', 'Lz']\n"
     ]
    }
   ],
   "source": [
    "print(list(output.simulations[0].thermo.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For backwards compatibility, the thermo can also be accessed as if the simulation was still a dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Step', 'Temp', 'Press', 'Lx', 'Ly', 'Lz']\n"
     ]
    }
   ],
   "source": [
    "print(list(output.simulations[0]['thermo'].keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot temperature vs. run step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd5gV1fn4P+/uwtKlK1VAsYANRESsKFiwxphYEjVqJInlq6ZKNHYTNbH8jLH3xFhjF7ErihRBERClSF1BqnRYlt3z+2Nm7s6dO+3evWV37/t5nn323nOnnJk5c97zlvMeMcagKIqiKF5KCl0BRVEUpX6iAkJRFEXxRQWEoiiK4osKCEVRFMUXFRCKoiiKL2WFrkBd6Nixo+nVq1ehq6EoitKgmDp16ipjTKeo7Rq0gOjVqxdTpkwpdDUURVEaFCKyKM52amJSFEVRfFEBoSiKoviiAkJRFEXxRQWEoiiK4osKCEVRFMWXnAkIEekhIh+IyNci8pWIXGaXtxeRd0Rkrv2/nV0uInK3iMwTkekiMjBXdVMURVGiyaUGsR34nTFmT2AIcLGI9AOuBN4zxvQF3rO/AxwH9LX/RgH35bBuiqIoSgQ5ExDGmGXGmM/tzxuAr4FuwMnAE/ZmTwCn2J9PBp40FhOBtiLSJRd1m7N8A3e8PZtVGytzcXhFUZRGQV58ECLSCxgATAJ2NMYsA0uIAJ3tzboBS1y7Vdhl3mONEpEpIjJl5cqVGdVn7vKN3P3+PNZs2pbR/oqiKMVAzgWEiLQC/gdcboxZH7apT1nKakbGmAeNMYOMMYM6dYqcKR6KrpWkKIoSTE4FhIg0wRIOTxljXrSLlzumI/v/Cru8Aujh2r07sDQ39crFURVFURoXuYxiEuAR4GtjzB2un14FzrU/nwu84io/x45mGgKsc0xRucKkKiiKoiiKTS6T9R0MnA3MEJFpdtmfgVuA50TkAmAx8BP7tzHASGAesBk4L1cVUwVCURQlmpwJCGPMJwT3xUf5bG+Ai3NVH0VRFCU9inomtTqpFUVRgilKAaFOakVRlGiKUkA4qAahKIoSTJEKCFUhFEVRoihSAWGhYa6KoijBFKWAUB+EoihKNEUpIBzUB6EoihJMUQoIVSAURVGiKUoBoSiKokRTlAJC1AmhKIoSSVEKCAf1QSiKogRTlAJC9QdFUZRoilJAKIqiKNEUtYDQiXKKoijBFKWAUB+1oihKNEUpIBzUSa0oihJMUQoI1SAURVGiKUoB4aAKhKIoSjBFKSBEA10VRVEiKUoB4WDUCaEoihJIcQoIVSAURVEiKU4BYaP6g6IoSjBFKSBUgVAURYkmZwJCRB4VkRUiMtNVtp+ITBSRaSIyRUQG2+UiIneLyDwRmS4iA3NVLzfqglAURQkmlxrE48CxnrLbgOuNMfsB19jfAY4D+tp/o4D7clgvTfetKIoSg5wJCGPMOGCNtxhoY3/eAVhqfz4ZeNJYTATaikiXXNVNURRFiaYsz+e7HHhLRP6BJZyG2uXdgCWu7SrssmXeA4jIKCwtg549e9axOmpjUhRFCSLfTurfAFcYY3oAVwCP2OV+Nh/f3tsY86AxZpAxZlCnTp0yqoQamBRFUaLJt4A4F3jR/vw8MNj+XAH0cG3XnVrzU85QJ7WiKEow+RYQS4HD7c9HAnPtz68C59jRTEOAdcaYFPNStlAftaIoSjQ580GIyNPAEUBHEakArgUuBP6fiJQBW7F9CcAYYCQwD9gMnJererlRBUJRFCWYnAkIY8yZAT/t77OtAS7OVV28aLI+RVGUaIpyJrWD+iAURVGCKUoBoT4IRVGUaIpSQDhoum9FUZRgilJAqAKhKIoSTVEKCAfVHxRFUYIpTgGhKoSiKEokxSkgFEVRlEiKWkCoj1pRFCWYohQQOlFOURQlmqIUEA5G3dSKoiiBFKWA0IlyiqIo0RSlgEigCoSiKEogRSkgVIFQFEWJpigFhIMqEIqiKMEUpYAQdUIoiqJEUpQCwkHnQSiKogRTlAJCFQhFUZRoilJAOOg8CEVRlGCKUkCoAqEoihJNUQoIRVEUJZqiFhDqpFYURQmmKAWEOqkVRVGiKUoB4aAKhKIoSjCRAkJEOovIj0TkYhE5X0QGi0ic/R4VkRUiMtNTfqmIzBaRr0TkNlf5aBGZZ/92TGaXExdVIRRFUaIoC/pBRIYBVwLtgS+AFUAz4BRgFxF5AbjdGLM+4BCPA/cAT3qOeTKwjzGmUkQ62+X9gDOA/kBX4F0R2c0YU123ywvHqBNCURQlkEABAYwELjTGLPb+ICJlwAnACOB/fjsbY8aJSC9P8W+AW4wxlfY2K+zyk4Fn7PIFIjIPGAxMiH8p8VEfhKIoSjRhpqKb/ISDzQBjzMvGGF/hEMJuwKEiMklEPhKRA+zybsAS13YVdlkKIjJKRKaIyJSVK1emefpkVH9QFEUJJkxAvCci7byFInI08GKG5ysD2gFDgD8Az4mVOc9vTO/bfxtjHjTGDDLGDOrUqVNGlVAFQlEUJZowAfEA8IGIJHphETnLLj8+w/NVAC8ai8lADdDRLu/h2q47sDTDc8RHVQhFUZRAAgWEMeYh4HbgfRHpIiKXA9cAw4wx0zM838vAkQAishvQFFgFvAqcISLlItIb6AtMzvAckWi6b0VRlGjCnNQYY/4tIluxopgWAwcbY1bHObCIPA0cAXQUkQrgWuBR4FE79HUbcK6xQom+EpHngFnAduDiXEcwgSbrUxRFCSMszHUGlhFGgBZAByyTkwDGGLNP2IGNMWcG/PTzgO1vBm6OU+m6ovqDoihKNGEaxAl5q4WiKIpS7wgTEItNxEwyEZGobeozDbfmiqIouScsiukDOy1GT3ehiDQVkSNF5Ang3NxWLzeoj1pRFCWaMA3iWOB84Gk7smgtVqqNUuBt4E5jzLTcVzF3qAahKIoSTKCAMMZsBe4F7hWRJljzFbYYY9bmq3K5QtRNrSiKEklomKuDMaYKWJbjuuQdVSAURVGCKcr1INQHoSiKEk1RCgiHBhyApSiKknNiCQgR2VlEhtufm4tI69xWS1EURSk0cVaGuxB4AStJH1iJ9F7OZaXyheoPiqIowcTRIC4GDgbWAxhj5gKdc1mpXKM+CEVRlGjiCIhKY8w254u9mpwOvhVFURo5cQTERyLyZ6C5iIwAngdey2218oP6qBVFUYKJIyD+BKwEZgC/AsYAV+eyUrlGJ8opiqJEEzpRTkRKgOnGmL2Ah/JTpXyiKoSiKEoQoRqEMaYG+NKbsK+ho05qRVGUaOKk2uiCteLbZGCTU2iMOSlntcoT6oNQFEUJJo6AuD7ntcgzqkEoiqJEEykgjDEf5aMihUAVCEVRlGAiBYSIbKC2L20KNAE2GWPa5LJiuUSjmBRFUaKJo0Ek5V0SkVOAwTmrUR5RH4SiKEowaWdzNca8DByZg7rkDfVBKIqiRBPHxHSq62sJMIhGYr43jeMyFEVRckIcDeJE198xwAbg5KidRORREVkhIjN9fvu9iBgR6Wh/FxG5W0Tmich0ERmY3mWkhyoQiqIo0cQJc33YGDPeXSAiBwMrIvZ7HLgHeNKzbw9gBLDYVXwc0Nf+OxC4z/6vKIqiFIg4GsQ/Y5YlYYwZB6zx+elO4I8km6lOBp40FhOBtiLSJUbd6oQ6qRVFUYIJ1CBE5CBgKNBJRH7r+qkNUJrJyUTkJOA7Y8yXkuwp7gYscX2vsMuW+RxjFDAKoGfPzDKAqJNaURQlmjANoinQCkuItHb9rQdOS/dEItICuAq4xu9nnzLf8b0x5kFjzCBjzKBOnTqlW43oEyiKoihAiAZhz6D+SEQeN8YsysK5dgF6YyX/A2vp0s9FZDCWxtDDtW13YGkWzhmAqhCKoihRxHFSbxaRvwP9gWZOoTEmrbkQxpgZuJYqFZGFwCBjzCoReRW4RESewXJOrzPGpJiXso1RJ4SiKEogcZzUTwHfYI3+rwcWAp9F7SQiTwMTgN1FpEJELgjZfAwwH5iHte7ERTHqlTHqg1AURYkmjgbRwRjziIhc5jI7RSbwM8acGfF7L9dnA1wcoy6KoihKnogjIKrs/8tE5Hgs30D33FUp96gCoSiKEk0cAXGTiOwA/A5r/kMb4Iqc1ipPqAtCURQlmKg1qUuBvsaY14F1wLC81CrHiDohFEVRIolak7oaaPBLiwahyfoURVGCiWNi+lRE7gGeJXlN6s9zVqsco/qDoihKNHEExFD7/w2uMkMDXxNCURRFCSfOinKNwu/ghzqpFUVRgomcKCciO4rIIyLypv29X8Skt3qP+qgVRVGiiTOT+nHgLaCr/X0OcHmuKpRPVINQFEUJJo6A6GiMeQ6oATDGbAeqc1qrHCPqplYURYkkjoDYJCIdsLNji8gQrDkRDR5VIBRFUYKJE8X0W+BVYBcRGQ90IoP1IOoT6oNQFEWJJk4U0+cicjiwO9YUgtnGmKqI3RoEmu5bURQlmEgBISLNsNJvH4JllflYRO43xmzNdeUURVGUwhHHxPQksAErUR/AmcC/gZ/kqlL5QvUHRVGUYOIIiN2NMfu6vn8gIl/mqkL5QH0QiqIo0cSJYvrCjlwCQEQOBMbnrkp5RFUIpY5U1xh6XfkGj41fUOiq5J3v1m5hesXaQldDySFxBMSBWAn7FtrrSE8ADheRGSIyPae1yxGa7lvJFtu21wBwy5vfFLgm+efgW97npHsax1hR8SeOienYnNdCURo4NRoRpzRC4oS5LhKRdkAP9/YNOd23g64HodSValswVNdoW1IaH3HCXG8EfgF8S63VvkGn+1YDk5ItHM1B5YPSGIljYvopsIsxZluuK5Nv1CqQPb5duZHHxy/k+pP6U1JSPCK4poFJhrWbt1FWWkKr8jivvlLsxHFSzwTa5roi+UR91NnnV/+eyr8nLmLeyo2FrkpeaWDygf1ueIfDbvug0NVQGghxBMTfsEJd3xKRV52/qJ1E5FERWSEiM11lfxeRb0Rkuoi8JCJtXb+NFpF5IjJbRI7J7HLSo4G921ln/dYqKn7YnJVjOaaWIlIegIbpe1izqdEZA5QcEUfPfAK4FZiBnfI7Jo8D92DNxHZ4BxhtjNkuIrcCo4E/iUg/4AygP9a6E++KyG7GmJykFdd03xan3DOe+as2sfCW4+t+sEQ/WVz3tiFFLzVEYaYUljgCYpUx5u50D2yMGScivTxlb7u+TqQ2K+zJwDPGmEpggYjMAwZjzbnIGQ3o3c4J81dtytqxnFtZbBpEQxIQ2dQcrn55RtaOpdRf4piYporI30TkIBEZ6Pxl4dznA2/an7sBS1y/VdhlKYjIKBGZIiJTVq5cmdGJvT6I8fNWsXJDZUbHUiycjrLYJiE2pFF55fbsKeT/mbg4a8dS6i9xNIgB9v8hrrI6hbmKyFXAduApp8hnM983zxjzIPAgwKBBg+r0dt7xzmxKBK58cQY927dg3B+H1eVwkVRV17Bx63batWya0/OkizGmzh27M5BuaOLBGIMxZBx51YAUiKwJs6rqdCzNSkMmUoMwxgzz+auLcDgXOAH4maldkKECayKeQ3dgaabniKyD/X/Vxm1c+aKlKi9ekx1nbRhXPDuNATe+k/PzpMu2LLzwzqTDhqZA3PnuXPr8eQzzVmzgyyXp5xVqSBrE9pC6PjZ+AQtimhy92rauq5LK9uoa5jeCiL5IASEiO4rIIyLypv29n4hckMnJRORY4E/AScYYd4/8KnCGiJSLSG+gLzA5k3PUZ16fvgyofy+Uk0+oLgRd0sJVm9iwtf6uL/XUxEUADL9jHCf/K/28Qn4+iE2V29m8bXud65ZtgoTZ1qpqrn9tFmc8GM/l5722KCH5wtQKPl/8Q7xKNhL+OuYbjrz9I5au3VLnY01bsja28M42cXwQjwNvYUUXAcwBLo/aSUSexnIy7y4iFbZQuQdoDbwjItNE5H4AY8xXwHPALGAscHGuIpisyuXsyLGoZ/KBygwExMMfz+e6V79KfHeuyXttR/zjw6SO9+UvvmPR6sI0djdL1mzmmDvHsTrAcbtuSxWPjV8QKcz9BET/a99ir2vfyko9s0mQacjRIDdVxnvltm1PvuYoJer3z3/Jqfd+GuvYjYUJ81cD2QkMOOVf4xn2jw/rfJxMCPRBiEiZMWY70NEY85yIjAawQ1QjW5Ix5kyf4kdCtr8ZuDlGnRs8NcZQUmgp5SITDeKmN74G4LqT+gO1WpFfhzl/Za1AuPzZabRuVsaM6/Iy1SWQRz5ZwOzlG1LKa2oMJSXC6BenM2bG9+zTfQf237l94HGCOsf6aHkKGulXVlnPP64bZntNcnsJi+QqVn+FScwLqj/veSaEaRCOiWeTiHTAdhrba0Osy3XFGjPZ7Du+WPwDs79P7ejSISsmJvt/WMfovDQbthbe/BKkGWyyzSerNlojv6rq8KcVZl75f+/O5b2vl3PgX9/l029XZVjT7BHkg3Cim+I66r33JOwe/FDkk/IauHwIFRDOpf0Wy0ewi4iMx5r4dmmuK5ZLCj1RLpux8z+691OOuWtcnY6RiYnJS62JKfjaGoJDd2OlJSCcHEtRI8CwZ3nnu3O44IkpLF9fya1jZ2evkhnivv9jZ37PV0utcZ4zQCiN2Ztt92gF1QH34A/Pf8mP7681LZ1673g2VRZ+cJAPElF9Ibd05nfr6HXlG/XamR0mIDqJyG+BI4CXgNuw5i08BAzPfdXqNzMqrIe7eHX60U/1zQeRrgbx2cI1KWVOR/nh7JV8G9Dgw6Jo8k1QTRztxun0Sl1vyHtfL09x0NbEvHX1YSC53TXy//V/pnL83Z8AtQOETDWIoISFz0+tYMmaWift54vXMmnB6rTq3FBJzAsKefIvffEdAO99vSIvdcqEMAFRCrTCciq3xPJXlAIt7LIGSzbUvuenWvP63v9medr7xhEQr0z7jrMfmZT2sf2oqq4JDeFMN8z1J/enRrs4l3TzmK8TDkmvNlFfNIjv1m4JnBi5tcoyt3g1iLnLN3DBE1O4+qWZSdsHjZ691AdTQ9D9T1eDqPJIxXSe65ZtxeWTiHNL60PbCCJsotwyY8wNeatJA8PpODLp89xmidUbK/nn+/P488g9aVpWK68ve2aata3tNI3L1EVrmPDtai45sm+i7O9vzebBcfN56/LD2H2nVNmeDZOX+xBbtlWnlIF/Z7plWzXNmpQkJuqNnfk9TUqFo/bcsc51CuLgW94P/M0ZTddqEFa91m6xQnW982Xi3ru6OivfmbWc0hI4co/M74vXuezgXHNpzHa23euDSKP9OAK4sRPnjtQ3S4IfcXwQjY5sXJjzvmfSubr3ueXNb3j804W8OXOZ77Yb04yn//F9E/jH23OSyhxb84oNW333yc7IvvYYLctLreN6NQhPx7Kxcjt7XjOWW978JqFt/Po/U7ngiSlZqE9mOB2Yo1Q5HbvTKXqFddz1IOqao+rCJ6dw/uN1uy9RGkRJnKB3Un0Qcc1sAFuzmO6jEIyd+T0H/e29SLOs057ri9acKWFN4qi81aIB4qjjmYwCnF0eGjef976x7I+btyW/OGV2j7Juc90nmZXab35QY3U6uTemLwv0H0Thd2jv+bw+iGmLLbPXA+Pm8/DHCzI6bxxe/uI7+l8zNlbIZWVVDUvXbkncE+f5OkK9zCsgYj7/qMCIlRsqWVKH2fzvfb2c6RXBZsR/vjeXj+f6R1I5UUxxTUxek6QzEFi5oTIy39PWqsxNTHOWb+Co2z/MyjuRKde+OpNl67ayamN47rbaqL7gBtIQljwOFBDGmFRPZCMhTt6hx8YvYLer3wz83RlJZqJBGPsduXnM14mJNFs8AqJFU2sUvj7DWchu+3+pfbnPfraET+eldhJvz7L8KBf/93OOuv2jlN/jmAX87oO3zPt99abal+ytr76PPEem3PD6LDZtq2bhqk2RQQWTF65h6C3vJ+ZIOHVOdKIeARF7hBjR5A64+V0OrcNCPhc8MYWT7gmeCX77O3N4/NOFvr+l66T2mpgcYXrAze9y0X/Cl6rfWlXNQ+PmM8Un0CGKO96ew7crNzG+gCHDZRGDLYfEwCJEHtZGOtVfY01MpbI4uOvdWtPM9a/NYtv2msCGUGtiSv88NcakjLS2VHkFhOUeWr8ls7BAd72cTu3Nmd9z1sOpju+gjgNg0epN7PGXsbwwtSL0fH5yMkqDcHc0pSWSNML/65ivA9NV3PzGLE7458eh9XHjvH4j7hzHYX8P74TneibPOaNjR8Pz+hLipk2pr13A7le/yUVPWZ2699qmV6zlwXHfpuzj9WVU15jEs3M04iA2bN3OzWO+5jSfQIconDkqzuApX1T8sJkV6y3zbJk92oobkRdnAFlf2wYUqYAIeiB3vTs3pVMLMkvUOqkz80GsWJ+sogaN0oMci1G46+198f89YSG9rnwj1nEck9Pr08NzJ/p1lO6qr1i/lVWeyCF3HctKhY2uCXQPjpvPQ+P8zU4PfbyAmd+tj6y7QzoDtB885gvjERApGkTM5z9pwRpOySDX0+qNlXy9LP61pot7DozXfHbSPeP565hvksq+WPwDb32VHLlXY0zKACeI+z+qFTgn/2s836/z94v54ZjInPW0X5++lF5XvpHovB1emfZdRokXvcxaup7VGys55NYPGPzX94Dae+T1w3hJ+CBc7WPszGUZaU6FpCgFRBgbPRN53PbWDVurEo5e513KZNF6AynJt7wmJkfwxD28V7C5v3s7tb+88hVe3Nv3uvKNxH1IV6VOOqarcPBf30tJhlflOub4eatTJuzFEY6H3fYBN7w2K2Kr+BJimqdjcR5/ZZWjQSRv770tYRqF99hxOPrOcRz3/+JrS3UhTqTVj+79lPc9WkKNMUntN25yxi+XrE3MBUiHMntyyn/sRIvzViT7zS57ZlpSW1uzaRvXvfpV2vN9Rt79Mcd67r3zPng1iOoawy8em5yYI+T87O4ffv2fzznt/gls217DOY9OZuZ39T8hRVEKiLD3wBnJO9tUuRrVkbd/xOCbrZGE8zJl4maqMYYZnsaR0tE45TUGYwxjZy4LHbV4TVZe800UXk3pux+sCU6OSh3l4PW7D1FCpcrzwno7ljgd1uI1m3l0fLiDuy4mXkdQOx2Cc0lbtlVz7qOT+dbTOWV7MmBQMsEgLnj8Mz4JcEZHETeKycvazVVJQRarNm6Lb3oLeTZD//Ye17wyM6XcObbTxstKwyt+0xuzePzThaF+riVrNtPryjeYvCB5hO+dL1Oa0CCSr2/5+q18OHsl//f0F0nlfs1hzvINjJuzkimLrAy3IlZCw9Pui5fQcMyMZYx+MT8r+hWlgAgjISDs7+5Zo05jGXTTu0yabzWkSQtWc+GTU9IKZ5vw7Wpe/DzZpl8iyaNPd5jc27OW8+v/fM69H6bag2vrHTz6zkRAOBEWTeyXr3J7DYNuepfnpyxJ2q6mxvKn+JqYQjqJHzZt44bXk0f+Sz3mhigB8cgn8SKf6mLjrXE9B4D3v1nBzW/MYvy8VXw0ZyU3vpF8Dd6OI4iXvqjgm++zbzp675sV/Po/UzPaNyiKKaqzn7xwTZK/yBgTmcPKwWvWcrN03VaenLAopdx51RxhHNS+e135Bi9MrXAlIww+1yRbMDwzOXilvKcmLWKhnYl4W7W/SW3Zuq386YXpiffHr1/wVuP612bxwtSKhMBw42ehuOipz3k6pJ7ZRAWEhy1Vyc5IxzG22hXWtmpjJYvWWA1l/LzVvDNreVppfS97Zhrfrkw2MT38yQJ6jx6T+O60i2pjEgnPKn4IjsDx+jDcI9k44Yvejs3pE5wGvrmymlUbK/nDC9OTtvvDC9PZ/eqxvhpEmIB4alLqi+/NnR8l1258Pcq0ZFEnDcKWm+77+dDHCxKjbe8lemcZe3E62yue/ZJj7/I3HW2q3M4Tny7MOIbePeES0jCDirBlW3XKYCGqHhu3bk9qf/+euCina3UnNIia6Al+t789O3E9Yds5SkiYT+mql2YmNCWvOdSdGfjZKUsSc36MMazfWsWTExYmfg8TVO7MDE9NWsStb30TuG0+KEoBERaT7thSnWfoqH3/nZQssb0vjdtXMeHb1Rn5JiA1bbZ7JnWYlSdMQMQJX0zRIJxRmt3Qg9Jx/M/WhLzzOKz6Bt8Dv9C+ZV4BUdfZZc656qBDJExMnusPCk2M0iDimKAeG7+Aa1/9ihemLonc1o+mdm+3aPUmLn36CzbHdCDX1Bj2vGZsSoqXqDrf++G3/Pi+2qikx8YvZG3MuQqbKquTBl+x6ulpm69OWxoY5FFjTKIdlpVIiq/Pwem0nW2jtKaqasPTkxcnwqbPe+yzpN/X2TPvq43hmpdnco3L7xc2YDn/8SmMm7MSsATSAx/ND61HrilKARHGloQPwnqKP2yu4qqXZqSMkL0vjdPwPvhmBWc+NJEHxmX2YJ0GWuOyeZfGiJjympi2bKtm3ooNVFXXxHLOVXmup8YzSgt6AcP68HQDsJyXyiFb4eHZ9EEkjumz7X8mLopM6x3HBOVE6UxbkpkT09EgbnhtFq99uZQ735kTsYeFM0iYOD/ZDu/cA2+0UBj/fH9urO3ufHcO+9/0bug2v3gseXFJ7zN5dPyChDbpHZTUmNrt5qzYwJ7XjGXMDCtrwaT5q7ltrDVCL/XMa4oSiuu3VDH6xRmc9fBE3983bXNm5BtWePwYUabTHzbXnxTpxSkgYjip3Tw1aTFfeJZM9DbE4XdYE8yW2GagW8dmpho6h02YeIxJNN7qGsNXS9f52q69Turhd3zE8DvGceI/P+HVL6OX9w5ygDsdWlBK8Jblwem8wtR1v+gV7zoR6eYvmrdiA72ufCMlOiSbPggHP+3o6pdncsl/v0gpdxOVGLGmxtDcjvOvjpCws7+3ZhZ7U6g4AqJ1M+vZxPXVBD1j51pH3Bk/rfxTk9KzkYdpmx/OXpn03U+rq7CDKryacE1NrQbhhL5+ONuKwjr9wYkJv16pR4OIEuROJ17xw5bQyazGpA7sohRjEYmV0SAfSxcXp4AIwck26R11ezsrvwb9wEffxnbOBVETZmIyhuPv/sTXdh2UwuCbmIsJBdmdHQ0iUHvxKf5hc6z1pPUAACAASURBVBUj7vgo9KX3hiYCzPLE+7tfJMeWG7Q27/frtibCL73Cpy4zVf18EJB5yoiq6pqkF3vZui30varW91RVU5M4V5Rp7JY3v+bblZuY8G1yCm3HxLRD8yZp1S3I1OM8R6+Gl00G3PA2N78Rz6eUMH+6nonzySuAa4xhua351CYlTO32Sj1m3KgQa7fPcf8b3wncrrrGpGjSUdpJqQgLY6xBnY/0+UUpIML6i6AJP95Oxu/h/O3Nb3jZ0zk9OO5bBt/8bmyfRMLE5HISJ0xMrmN8vWx9UtqIuiZB+2hOsmnEuT5H4AXVf0PAAjBzV2xMe8F2rzBzC+UnPl3IPte9Hbg278dzV9Lcnn2+fP1Wel35Bu/OSj8Vu5cgDSJurL+X7dUm6Vjvfr0iaVAx+/sNTLdNS89O8fdBbNlWzZSFawLbatOyEv7y8kye8IkACmN9wEp/1TWG976u+72MOvdDIfm43ELVeSbuMYvzuzd0usZYbRFqrQNPT16c8p7WCghr/0gNwiUgwgaFNcakDK6mR5gOS0viLdUaN2KuLoSl+y5K4trag0bHbvPPja/PSqj3UdEtDm999T2nDuxe2zEZk4iwWOKKYvJOnqrrsqHeiKDt1TU88smChPkik8HKOY9Ojt4oBPd79XZEZ9+5TbPETO3pFdYLeO+H8xjer25pw2vNGck3wDuhMi5V1TVJgwvvcwvLpwTWTGTvbGbvsp5lJcK/J6YnHNx4B1DVxhQ0wy4kt4V3Zy1PTFjz/p6ymJFrR7cJ7fJnp7n2dWvpVlnU6HxNTCf8qH+nhhyv2hTulBeRwH7IPfO8qqaG5uQ27UhRCogwxf2N6ct8Vek4JiZIbqBu269fVEdpiaQc57fPfcmpA7snGvzUhT/QrIn1MoSll8j2aKKq2iQJjUKkLd5eY/hyyVr6dW3DxPnhK5GVl5XQrIn1sjjLWpaIsHTtlrQXRHJTq8l5J/VlJiCenryYrm2bJ76nK9hveTPVt+VND1JX574xJDnbr3B1poXC3dE/MWERT0xYRPd2tffRYFi1sTIlf5dbsHyx2H8mu7tpz/5+PRu2VkWamOqy1nbUu/qXl2fyk0HdfX+b7ErV4U2fnwuKUkCEMWH+aibMX81ObZrxvStqI9NZpg5+6YFbNi31VevPeHBColMLMjN4yTRnUxAbK5M7nVzGtQexaPUmbh37DWcd2DNSg6mpqVXl3aP7oSGLA8Wh1heTHQ3CO9mxrpof5GaNhbMeqg11HT8vP8uENi0tYdbS9ezZpXWKSXeAj53fvUnVdsOgm95l3x5tk7aJY6r59NtVVNip1pevr+TnD0/inrMGhu6TaZZliM7jtGJDJf/6IHVSrDGG1q6gkLhWibpQpD6I6CGWd3H1uqbk9ZtI52Rs9eINM4xDNjoaN14tqhACwplb4Z2D4ke1MQkh6XQKmY7y3QT5INZnyWEbNCM3HT6IyKDaUNhWXcPIuz/mlWmpUXd+z9LtxK+0n7k3SV9QZJabsx+ZnJSf7MuKdbwTYdIMmk8Rh88DNJkoXpm2lNdcSTPzodXnTECIyKMiskJEZrrK2ovIOyIy1/7fzi4XEblbROaJyHQRCRffeaDSI+XjLqYSxPL1qRqEYzrKBtmOaPCmGS/Ewlhr0+iEq2sMVduTHfyzl8eL4ArDL2IG4keHRVGXBXQc5izPbJGn+srlz07j+LvTS1CYjeytbrxpYLz4TQyNyyc+a7LE4fJnp/Hi57XO9Xw4qXOpQTwOHOspuxJ4zxjTF3jP/g5wHNDX/hsF3JfDesWKi08Nc63bOf1C+BybeTaIUlvTxRulUwgNwplRGofx81blROWurjHMWro+xeacrRTcb8/K/kJJdZk5Xl/4amn0/S3kOjtx05vnkjjms7qSMwFhjBkHeG0lJwNP2J+fAE5xlT9pLCYCbUWkS67qlgl17R69jkTIroDYluXRhNc3Ut8XWH/o4wWR4YOZMHnBGkbe/XHKQjjZ0tiWrEkvFDguTUrrj5Bo2yK9+RhxqcsSrXUlndxruSIfJqZ8O6l3NMYsAzDGLBORznZ5N8Dtja2wy5Z5DyAio7C0DHr27Jnb2rrIhbTOponpxtdnZXXZzmzZ2PPJ0nXZ72zHzExpgg0Cy2dWP6T6jq2bxc7NlA6FMHvWJ+o6KTcO9cVJ7Tfc8b16Y8yDxphBxphBnTp1yuxkGQyusu0EBmieRQ0CSMllH4ddO7fyLffmj2kIfJzhOgjd2jbnkz8N4/DdUttTXZyRhWLRmk05aa+ZsuMOzQpdhXpPnJT8Xhq0kzqA5Y7pyP7v6O0VQA/Xdt2B6ARCeSRONES6ZNPEFIdP/jSMx35xQOL7kD7tOWnfrv7bZuhIa4iIQPd2LXz9LGGmpO7tmvPiRUMTqS3qC34BEdki3fQdADu1Kc9BTfJH2JoV2SKTUzTGMNdXgXPtz+cCr7jKz7GjmYYA6xxTVC7IxImXixFZvgVE93YtGLZH58T3Z0YdlLEt9cZT9uK/vzwwW1WrF6Q7IuvQqpyBPdvRv1ubHNWo/uEkAEyHHds0bA0ik2tOl3QTU0Jmyx2nSy7DXJ8GJgC7i0iFiFwA3AKMEJG5wAj7O8AYYD4wD3gIuChX9cqUsNm43gVa4pJvAeFHphkhm5YKQ3ftGLrN4+cdwN9P2yej4xeCdB3PzqQl74u6/87tslan+oY3LLeba1Z4EPt5Jq81NMIyFmeLTAREgzYxGWPONMZ0McY0McZ0N8Y8YoxZbYw5yhjT1/6/xt7WGGMuNsbsYozZ2xiT08QvcZ9Fx1a1qnHQNH2AFk0z6+iz6aTOFLfpbPKfj4q9nzNxcM5Nx/n+Xl5WwhG7d6Y8TSF4aN9woZMLnPaQ7ojMWbfB7Sy86/T9uPHkvbJWt/rGz4ckB4b8/bR9+OOxuyeVuU1u/73wwJxFMWWLs4fsHPp7kzyYEDMxMYWl088Whe+h6jGXDNsl1nbNypI7wd13bO27XafWybbYTDWPbOIWEJ3TMAU47TnoGpymm+4Ewz4dW6a1fTbo3dFy1MfVIBybdCvb9PCLg3slfispkXoVYpptzj+kd9L31s2acFjfZOd+uatNlIqkJNYrBBcP2yVQm/U6iI/fOznC3vs8d4zpU7nttH2Yft3RzL3ZfxDlJjMTU9q7pE3hn1w9Jq4JqLlHg/A+6+ZNSrnvZwO56/T9ksrz4fyKwllo6OBdOwDJJoMBPYNNA5GDF/v3IX3ap1WfZhlqY5nw9IVDeODs/fnnGQOA+JMBHUHiaBA/HVQbX7Frp1Z5GXHmgs+uGs4x/YOz3546oBstPelhykolpYMtd2nGpSVCWT0QmMN278xPBvWI3hA4Y3Dydl4B16FlPAHRqVU5bZo1idUegtLmh5GPyasNsyXXkbjCOm7oWSuPjdKbt2ngzm05bu8uKa5xv4VL8o3jfD/3oF4AjPvjMDq2agpAl5DwRG8qEgdHAzC2hOjQKr0Ilp3bR2sQ3hFeprQqL+OY/juxg20C8dp0Tx3QLel7uUdb8nNe9uvahi5t03PKnjqgG/f/PLfZZXbu0CJym06ty0Of10n7dU0xhZSVSMpAx93+LY0q83aerdnS3vXNW7oGIt5zeH0OXg3Cz/bvd38zCV1NBzUx5YjyslJuOLk//70wPAonrtrnOCUdE1LQXl7BUSgfxE5tmiXO7ZiYHF9BaYkkHJGdWwd3dH5RXTeeshd3nWFpSe536M8j9wj007x2ySHJddshWqBka9Ki9wX2Pu/fHJFsYrx42K788pDenLCPJaC8AwOH8rJS9kojsunE/brSr8sOsbfPhFGH9QGiZ1iHabUtmpaltOGy0hLKQgRAqdRNQGSCX1vzPtsnzh8ce3/v9fl1zH73LdeaU4OOYqrvnHNQL3buED5ajTsCaFleyvy/juRiV4fyv98clPjspEv2Hq51HqIj/Pj4T8OYcd0xAJxn28/36lrboTn1DAtP9BMQp+zXNREn746OGnXYLnx1/TG+x3F3pB/94YhY9upspbnwvsD3nDWAX7ps7B1blfPXH+2dMLvt1KYZV5/QLzHBMSy65flfDeX93x0eqx4tmpSmjGK92oqbNy87NO3oMCfgwu1QP3VgNx44e/+k7cI68zbNU6/XT4NwU+rz++jj9mD4np0D9kgmk0GyX5Zkry9sUK/2XH38nkBq2Lt3AmscDcKv3eba99Kgo5gaAlH9f0mJ+DpNU9XsEkpKJDHSEIH9d26fEgLoHX21Comv7tS6nMlXHcWNJ/cPr2QGNCktSXQER+6xIwtvOT7JtPD0qCH8dsRuofWr9FmDoElpSUKoepuu+9o//P0RvuU7d2gZyy/jaBB19eF4BwA7d2jJ1Sf0S3zfoXkTzjqwJy9dNJSfDurOSft1TTp/2Ez45k1LA9O5e/Fz9O/VLVij2LNLGwb1SvbtRPl63IvrONzx0/0Y4AlBDbqlpw/q4Rt84eeD8P7uFTq/OnwXHj73gIA96k6JwIn7duUo15yfdPrqrp731lt/v7VX/LSFXJuY8pFqpMgFRPgDLBF464rDUsq9pgXnBfeONFJNGMnH6dvZP9oJLIdh59bNYncy2aR/1x34v6P6hnbAfhqENZq07kXYyK9XWKRSjHdq1GF9aN2sjBExlhP1hmC6CYqwcsIeHbt15zbNuO20fRNBC86cmKgotLhRamUlJal+LOd/wP3wPpu2zZsmffeaUBzzZ+tmZfzq8D5Mu2YEkGqbD+LXR+ziuyZKWUlJSufo/tarQ8u8O6mrqmv455kDOHVg7apsfp2100a9l9WktIRhu3dKRCt577Vf9JCf5pXtaLZeHj+HOqlzjFtAnOiTciLIftrGk27AaQheldKrvrtfsEl/Poo9u0Tbqb0RUpmQqaMvbATkJyBKS8JHk3Hwtnk/IbB3tx2Ycd0xKSO9Pxyze4rfoEPL5I7TTVBdbzxlLxbecnzgftvsdSeibOtxO4iSEmjXsinv/e5wfjtiN6D2mQUJMa/w8c41KCuRpGvo1Kqc0/bvzuPnDWb0cXvStkVT3+MH9Tk7t/d3cjcp9Q9j/e+FB3LjKXvRrEkpTbJgavETtj8d1J02tpZ71+n7cfnwvkDtyNr9fMPCrd2/fPzHYQA8dt5gLjnSOp73OfuZdvyedTrvwvA9owc7e+yU3F+oiSnHuJ+f3wQtp0P3jiBaN/MKCOs2OiMlpy16BYC7je7YplmsBpQNATH16hFMHB1/EpxDqAbh4ygWCbdHx8Hb6N2TrvbpbpldHFOet9O4eNiuKQIhbBSf6cjWufYwP4H1e7xn18eeh7FLp1a1bcjutrwj/N/ZAsR7n72DFm/bEhH+8ZN9U2Z5x9UggrYLGhQM3aVjQhNz3+dLhu0a63xeygOEsfOODt2lA71sn6Lj/3LfIz/tx4m0E7H8S5BsXnJ29woIEbjvZ8lRZ849uPbEWhNluxbBg5Pbf7Jv0nfHRHjm4J7MvP4YendsmWQiu+On+3LH6cn7qAaRY5JMTD732nno3hfAG97oOGadhuR2er12ySGJUYlzPvdpX73k4JTzPnzOoMTnbGR8bd+yKTtlkFEzLDrl2P47+ZaXhnS6Q3fpkPg89vJDueXUvVO2CQvde+Ds/Rl7+aEJc4xfB+19mZuWBt+/TFcJdOYKBGXCTZy7rIRp14ygf9dgTbFEAgYBARrE8XYEVRPPtXsFRlxBHTRI+dOxe/DZVcMj97f8WV5hlLyNW0D8/phak98p+/knivQ9j8+zNqa2Ey8pkURkntOE3G0x3MQkPP/rg7j9J/v6ah1++x63dxc+/uMw+nVpw9H9dky0uz6datuEd2KsG2fg0rN9CyZfdVSiLi2altKqvIwPfn8E151U6388dWD3FHNzPgREYcJo6gluAWF8JIR7BOGecdzS80I7zmjnpXS/IHt3r3U2OsdzN7d9uic7Cc8c3JPhLrOKIyBaNC1NLHPYZYdmLFu3NezSsoJfJ/O3U/emU6vywDxMYR3TE+cPTjh499ipTYrKDOG5odq1aEqXHWpHeH7agVcrCNMgMjWHnTW4Jz8a0C2Wf6hti6Y89csDmf39Bk5/cCJgzbAtLyvhsmempYxsE51WQB0T/i6P2aa9R3OKe21BQrK0JLyDc59HIvqpIBPT33+yLyfs05VfPhmdWScoY677/jmh2k51ynw6ezfuavdo34IeHjOa0z8E3coe7Vsw5rJDATj7kUkpv4eZIB2h2qq8jM6tm9VqM+7zRzzDPCwoV+QCwvX8/Pol5wEdv3cXnp1Su55RC3sEe8YBPRjQsy2De1vqYZRN2omvDjIdtm5WxrlDk/PCOB1C2+ZNEgJi7GWH8fsXvmTt5m18tvCH0HPWBb9O5szB4Ys0hXVM7ugpN/f+bCBtmjmT1YKP7dUY/DqNFA0izMSUoW1cRNIKHmjbommijQCcNrA7qzZZKbm9d+uUAd14/NOFnDu0F5MWrKF9y6ZsdM2ydWYyu0ftj593QEoa7rjX5t3MaZpxMx6XlQg1JiLYI6BNNCktoU+neKlVgp6jc2hjalPeJPw3rvM61/nSRUNTBF9Q7Z3j+Jmnwnj90kOo+CF88Spnkqz30O7vURquzoPIMckaRPDvN/2oNvlaaYnQ2W5grcrLOP2AnrW+ioT92J82zcKTls247piUUbXTEezrCkfcoUUTHjpnUOhEtmyQiT8hk0535N5dOMT2ATl+Bj/zg/dF9TMxeescppGEmcOyjXd2camPuREsbfSzq4Zz3F478ZcT+iUmHoKVGqSdrSm4O78jdu+cIrCcx3Dewb04/+DegfWKY2bbOyTkViReYMJFR+zCC78+KKU8zE4PtZNJgwW9dW5jTGJbJ9rMfW1OHQf0bEf3dpam4LzHQRkD0smPtIttWmrXogl7dduBY/fyN8E6OH2Fc12+A9SI06sPIsckCQg/H4T9u3tUOnH0UXS1zRzeXRKdY0DD8joS47Br59Y8ePb+3OYzMSrXcdZhPogg6lqlHds0Y+Etx3OkK6rjdyN2843BP27vLvS2Q2adW+59qbu2bc6vDuuTWC3OiRKCzH0Qftz2433SSpfhPLugkbqIcMEhvROdGMBBLh+OV1h6Z/869+HaE/tzjctxGlQPB2/o58TRR/Hsr4aEXEk8/njsHilzNyD1nbhqpDV5rWlZCTee3J8De1vX7Kd5GuDmH+1F93bNadeyaUIwOILCPcD2e9Y/GmBNFDzHTjPjd3yIFXnN6JF78OT5g1NMxkE4AxlHC06cyzOQCCMfqTaK2sSU7KP28UH49I+dWpcn9vM+n6iwxqiolyCO7r+T70g418n+4h7/9UsP4bu1lkrtNPCRe4ePoKJwRnXrt1Zx6VF9ufSovinbdGxVzge/P4KnJy/mALvz6de1DXt1a8PM79YDVgc4euSe3PDaLD6as5IWTa1Zy8ZkV8D+9IB4ieAcEi9/RBXiz+b35gOLV48o80lYcMPrlx7iW57OglylJcKNJ/dncO8O7L5Ta1ZvrOTmMV/TqVU5Zx/Ui7dnLQes9UccbjplL65+eSZnHNCDQb3ac4wdMOFcihM95n5n/K5TRBL7+rHJNu3FWQ+ivKyUw3yWrHV4+4rDEGDEneOA2mwAXg3CXcv6YGIqagHhfvn8VTz/B+SUe1U8pxEG2zStXzIRFH4NPNcaRNzj79Vth6SZv9OuGVHnRVb279mOUYf14eQYkS5uv8huO7bm9UsPpdeVbwC1L5nbVu08tkJm042rvcTdLtM1SepC2GzvdDjbNYJv37Iplx65KyPthIw/HdSDj+euol/XHfiyYh3/d+Su/HzIzvzcZw2HFk2sNufMDk/SIDJ41j9stlZb9M4xyaTV7ObRgJ15RI4w62bPdHc7yv36H2dwA/mZSV3UAiLKxhjUqGo1CO8TMkm/+/HiRUOztgRjUP3++8sDWbRmM6NfnFGn4zsdaImk1xjbRtiV41BSIvzZNjfUhYQpx75VNcZw22n7cNc7c2LPAcgFCQdoxHZxO7ZMtVMvlx65K2s2VXJGRDBCrhARfnd0bSjsift25cR9u/KvD+YB4WaVnh1a8NyvDkr4sdwDuEzMiWs3VwFWe+7XpQ29O7bkjRnZWQm5VkBYz+3EfbrQrkUTDnFFB/qFPwu15iidKJdj3O+ev5Pafz+nOJPHM7Bnu1jLNMbBiW658/R9kxYlGbprx8hoozg4Poj6sOBLppQkBITtzMQalX6awcTBbOKMHEcft0fodmECokRqnfnpRtoE0a5lU+46Y0Bgpto4ZNG1k+C0/bvTp2PLyHY9uHf7hC/CLUwyacK72PNc+ndtw5jLDmX0yPBnlQ6OTyWRAVqEQ/t2SnqOjvkpSDvUeRA5Juml8rnZgSYmJyFdHlS8MBwVuEPLcpqUlvDltUdn9eV0NIjSEoHU3HwNglKP2a/Qz8yh1JMKI2y7IOb/LXn/py8cwmXPfMGKDZVp+QEaAju2acb7riSPcejoWtgnExPTzw/syYG926eYh+rKoX07cnS/HbnmhH6RAu+Lv4xIeqfFZWPKhwZR1ALCjdvsU1oiVNeYQAHhlHoluJ+jKZc4ERCd7aRi3lj4uuK8VPVh5btMKU3RIOqJhIhJOh3bQbt04KWLD+b5KUvYbcfwWd5e0lm/wo9xfxjGxsrtjLz74zodJ5skT1JNvw2LSJJwcEKJB/RsF7RLJLNuOIYmpVb2Z+/yrX6080yAfPE3Q3luyhKemrRYfRD5xJ0UzhEQwT6IWnOFmz26tKFz6/KkdAK55NbT9mHszO8D18CuKwkNoh4sGZkpXh9EfdEg4pJux9atbXMuH75b9IYuJow+ss6Di54dWrB0bfjksEKSiYDw0r5lU1695ODQLMxR1DU787492rJHl9a2gFANIm+4zU1lJcI2ghvVgfas2OM8k2FalZcxOUb+mmzRuXWzwBjubFDrg2jAAsJjYspHaGA2yce9d6cvqQv1+c5mK+Iv7jyHXOK0aTUxFQinMbkdWx1aNuVAO+Ni3x1bx7Ifx2XM/x3quwhJoXEaYsN2Utv/A7S++k6uQ5lzQX2scQO8jYE4baLRCggRuQL4Jdb7OgM4D+gCPAO0Bz4HzjbGbCtE/RKmFZcGMfUvI3J2vn4h2T7rAw2xk3JICHtXmGuhOHy3Tny1dH1a+2QrOimbfJHDdyFX1Mf7mCkJE3djNDGJSDfg/4B+xpgtIvIccAYwErjTGPOMiNwPXADcl+/6Aa6lQxtPo6oLDfk2JMyEiZeqcHXxrvKWDqMO65PFmtQNr+PUIR8dlmJRWiKNOtVGGdBcRKqAFsAy4EjgLPv3J4DrKJSACFgHotBMHH0UW6ryF29qYkz8q+/U5jyyaIhdWDbNmcVGvlLj55tSkcaZ7tsY852I/ANYDGwB3gamAmuNMU5e4wqgm9/+IjIKGAXQs2duZntG5YEvFJks+lMXuuzQnD12as2fR+7JOY9Ozuu5s0VK1lQd5RYVr1xyMPNXbip0NbLOvT8byM4d/JeBzSaFMDG1A04GegNrgeeB43w29X2TjTEPAg8CDBo0KCdve33THDLl2P470aVt5kKlaVkJYy8/LIs1yh93nzmABz76NjGpsTZ/ViFr1bjZqU0zhvRpn3aYbS7p3LpZztPiF4LhPmu154JCmJiGAwuMMSsBRORFYCjQVkTKbC2iO7C0AHXDqpP1v6F3JvefvX/WjnXGAT04YvfgbJX1jZP27cpJ+9Ym+qs1MTXwh1qPKSst4ZlRqWs+KA2XQgiIxcAQEWmBZWI6CpgCfACchhXJdC7wSj4q89cf7c3uOyXPOq21Rmhn4nDLj1PXo2hInH5AD96c+b1vFlBFUfwphA9ikoi8gBXKuh34Astk9AbwjIjcZJc9ko/6nHVgqh+jocbMK8F0btMssX6woijxKEgUkzHmWuBaT/F8IPM4wCzgRDy0KM9/bn1FUZT6hs6kdvHapYewZM1mOrUu539Tv6NPx3gLqiuKojRGVEC46NiqnI6trMyolw1PXeJSURSlmGi4SXYURVGUnKICQlEURfFFBYSiKIriiwoIRVEUxRcVEIqiKIovKiAURVEUX1RAKIqiKL6ogFAURVF8kYackE5EVgKLMty9I7Aqi9VpCOg1Fwd6zcVBXa55Z2NMZHrmBi0g6oKITDHGDCp0PfKJXnNxoNdcHOTjmtXEpCiKoviiAkJRFEXxpZgFxIOFrkAB0GsuDvSai4OcX3PR+iAURVGUcIpZg1AURVFCUAGhKIqi+FKUAkJEjhWR2SIyT0SuLHR90kFEeojIByLytYh8JSKX2eXtReQdEZlr/29nl4uI3G1f63QRGeg61rn29nNF5FxX+f4iMsPe524Re5HuAiMipSLyhYi8bn/vLSKT7Po/KyJN7fJy+/s8+/dermOMtstni8gxrvJ61yZEpK2IvCAi39jP+6DG/pxF5Aq7Xc8UkadFpFlje84i8qiIrBCRma6ynD/XoHOEYowpqj+gFPgW6AM0Bb4E+hW6XmnUvwsw0P7cGpgD9ANuA660y68EbrU/jwTeBAQYAkyyy9tjrQPeHmhnf25n/zYZOMje503guEJft12v3wL/BV63vz8HnGF/vh/4jf35IuB++/MZwLP253728y4HetvtoLS+tgngCeCX9uemQNvG/JyBbsACoLnr+f6isT1n4DBgIDDTVZbz5xp0jtC6FvolKMDDOQh4y/V9NDC60PWqw/W8AowAZgNd7LIuwGz78wPAma7tZ9u/nwk84Cp/wC7rAnzjKk/aroDX2R14DzgSeN1u/KuAMu9zBd4CDrI/l9nbifdZO9vVxzYBtLE7S/GUN9rnjCUgltidXpn9nI9pjM8Z6EWygMj5cw06R9hfMZqYnEboUGGXNThslXoAMAnY0RizDMD+39neLOh6w8orfMoLzV3AH4Ea+3sHYK0xZrv93V3PxLXZv6+zt0/3XhSSPsBK4DHbrPawiLSkET9n0iy8qQAAA+VJREFUY8x3wD+AxcAyrOc2lcb9nB3y8VyDzhFIMQoIPztrg4v1FZFWwP+Ay40x68M29SkzGZQXDBE5AVhhjJnqLvbZ1ET81mCuGWtEPBC4zxgzANiEZRYIosFfs20TPxnLLNQVaAkc57NpY3rOURT0GotRQFQAPVzfuwNLC1SXjBCRJljC4SljzIt28XIR6WL/3gVYYZcHXW9YeXef8kJyMHCSiCwEnsEyM90FtBWRMnsbdz0T12b/vgOwhvTvRSGpACqMMZPs7y9gCYzG/JyHAwuMMSuNMVXAi8BQGvdzdsjHcw06RyDFKCA+A/rakRFNsZxbrxa4TrGxIxIeAb42xtzh+ulVwIlkOBfLN+GUn2NHQwwB1tnq5VvA0SLSzh65HY1ln10GbBCRIfa5znEdqyAYY0YbY7obY3phPa/3jTE/Az4ATrM3816zcy9Os7c3dvkZdvRLb6AvlkOv3rUJY8z3wBIR2d0uOgqYRSN+zlimpSEi0sKuk3PNjfY5u8jHcw06RzCFdEoV6g8rMmAOVkTDVYWuT5p1PwRLZZwOTLP/RmLZXt8D5tr/29vbC/Av+1pnAINcxzofmGf/necqHwTMtPe5B4+jtMDXfwS1UUx9sF78ecDzQLld3sz+Ps/+vY9r/6vs65qNK2qnPrYJYD9giv2sX8aKVmnUzxm4HvjGrte/sSKRGtVzBp7G8rFUYY34L8jHcw06R9ifptpQFEVRfClGE5OiKIoSAxUQiqIoii8qIBRFURRfVEAoiqIovqiAUBRFUXxRAaEoaSIiV9kZR6eLyDQROVBELheRFoWum6JkEw1zVZQ0EJGDgDuAI4wxlSLSESsz6KdYMeqrClpBRckiqkEoSnp0AVYZYyoBbIFwGlbuoA9E5AMAETlaRCaIyOci8rydOwsRWSgit4rIZPtv10JdiKJEoQJCUdLjbaCHiMwRkXtF5HBjzN1Y+W6GGWOG2VrF1cBwY8xArNnQv3UdY70xZjDWLNe78n0BihKXsuhNFEVxMMZsFJH9gUOBYcCzPiuTDcFatGa8vZhXU2CC6/enXf/vzG2NFSVzVEAoSpoYY6qBD4EPRWQGtQnQHAR4xxhzZtAhAj4rSr1CTUyKkgYisruI9HUV7QcsAjZgLQELMBE42PEv2NlJd3Ptc7rrv1uzUJR6hWoQipIerYB/ikhbYDtWJs1RWEs7vikiy2w/xC+Ap0Wk3N7vaqwsogDlIjIJa4AWpGUoSsHRMFdFySP2okcaDqs0CNTEpCiKoviiGoSiKIrii2oQiqIoii8qIBRFURRfVEAoiqIovqiAUBRFUXxRAaEoiqL48v8B+8ZXncpxMn4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Pull out data\n",
    "steps = output.simulations[0].thermo.Step\n",
    "temps = output.simulations[0].thermo.Temp\n",
    "\n",
    "# Plot\n",
    "plt.plot(steps, temps)\n",
    "plt.xlabel('Step')\n",
    "plt.ylabel('Temperature (K)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot pressure vs. run step (note unit conversions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEICAYAAABBBrPDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd5gUZfLHv7WRnBckL5IUUEBWBBERCYJ6YBZMmM9T7wx3ehjPM5yenvH0VMz+zjOcEQVEQAQPyTnDSs5Lzmyq3x/Ts8z29vR0eDvMbH2eZ5+d6el+uzq8b71Vb731EjNDEARBEOKRFrQAgiAIQrgRRSEIgiCYIopCEARBMEUUhSAIgmCKKApBEATBFFEUgiAIgimBKgoiGkREK4kon4hGGvzegogmE9F8IlpEROcHIacgCEJlhoKaR0FE6QBWARgAYBOA2QCGM/OymH1GAZjPzK8TUQcAY5k516zcBg0acG6u6S6CIAiCjrlz5+5k5hyj3zL8FiaG7gDymXkNABDRJwCGAlgWsw8DqKV9rg1gS6JCc3NzMWfOHMWiCoIgpDZEtD7eb0G6npoC2BjzfZO2LZbHAFxDRJsAjAXwe6OCiOhWIppDRHMKCgq8kFUQBKHSEqSiIINtej/YcADvM3MzAOcD+D8iqiAzM49i5jxmzsvJMbScBEEQBIcEqSg2AWge870ZKrqWbgLwGQAw83QAVQA08EU6QRAEAUCwimI2gLZE1IqIsgAMAzBat88GAP0AgIhORkRRiG9JEATBRwJTFMxcDOBOAOMBLAfwGTMvJaLHiWiIttsfAdxCRAsBfAzgepZ0t4IgCL4SZNQTmHksIoPUsdsejfm8DEAvv+USBEEQjiMzswVBEARTRFF4zPwNe7Bk876gxRAEQXBMoK6nysDF//oFALDumQsCliT1mbt+NzLS0tC5eZ2gRRGElEIUhZAyXPr6dACilAVBNeJ6EgRBEEwRRSEIgpDklJQyZq/b7Vn5oigEwUMKi0sxafn2oMUQUpzXJufj8jemY+aaXZ6UL4pCEDzk+QkrcdMHczD9V28qsCAAQP6OgwCArfuOelK+KApB8JBNu48AAHYdOhawJEIqk54WybFaUupN4gpRFILgIaTlSPaqAgsCAKRpL1qJRxmORFEIgodEe3qlkqJM8JB0rSUvFYtCEJKP9GhPrzRgQQTPOHisGDe8NwtvTV2D3JFjUHDAfzdjmevJow6JTLgTBA9JE4si5fl24RZMXlmAySsjKyDk7ziInJrZvsoQdT2JRSEISYimJzyrwELwRHvz8b77KYMMZguCRbbsPRK0CGV47RIQgifqXowSgJ6IGcz2qHxvihWE4DjzmR8xYVk4Jrl57RIQgicjvbxmIArOohDXkyDYwMt0BnY4HvUUsCCCZ6SFwKLw2nIVRSGkJIcLi4MWAUCMS0A0RcqSkaZXFP5rCq/fs0AVBRENIqKVRJRPRCPj7HMFES0joqVE9B+/ZRSSk8OFJUGLACDG9SRjFClLWigGsyP/U871RETpAF4DMBhABwDDiaiDbp+2AB4A0IuZOwK423dBhaTk8LFwKIqyCiyKImXRWxR2DIriklKMXbwV7PL9SE/hmdndAeQz8xpmLgTwCYChun1uAfAaM+8BAGbe4bOMQpJyKASup2venom3fl4LQCbcpTJ6i8KO6+mtn9fi9o/m4dtFW5XIkHIWBYCmADbGfN+kbYulHYB2RDSNiGYQ0SCjgojoViKaQ0RzCgoKPBJXSCaOhMD19L/8nWWfxaJIXdyMUWzbpyWNPOhuNncqWxRGd1N/lRkA2gI4B8BwAG8TUYUFkZl5FDPnMXNeTk6OckGF5KM4ZIPHEh4bXpgZa3cecny8fh6FE75btBVjFzuzKtbvOoQ9h4sAeGe5BqkoNgFoHvO9GYAtBvt8w8xFzLwWwEpEFEfombNuN3JHjglajKRj4+7DWL39gOtyAgg8MUUm3IWXd6etQ99//IRFm/Y6Ol7veuIK/d3EzF2/B7d/NM/R+fs89xPenRZxcXpluQapKGYDaEtErYgoC8AwAKN1+3wNoC8AEFEDRFxRa3yV0iEfz9qYeCcLHC4sxqezN7ge7EoGDhwtQu9nJ2PAi1ODFkU5YlGElznanJshr07D0Ff/57q8IKtqyoXHMnMxgDsBjAewHMBnzLyUiB4noiHabuMB7CKiZQAmA7iPmZNiqbBjxWp85E+OWY4/f7G4nL87FVmwcS9OeeyHoMXwDLEowktsL3zhpn22j9c/2lRUFIFmj2XmsQDG6rY9GvOZAdyr/SUFRSWluP/zRVilwH0CADv2Rwa5wjIvwCucmv3JghgU4cWtX1/varLjelKd7sMr15OkGVfM7HW78dX8zQpLjDz4kLnclZPq1yeup/Di2q0bIovCKySFh2JUT9+PvnRBJBpLZrbvO4q3pq4JzdiOpPAIL2574ZXhyYqiUIzq5jz6EgaRaMxXFCvCLfuO4qmxy7Fh92Gl5TpFxijCi1sdrnKM4tznf8KbU351fLxXzYQoCsV45XNcsnk/fi04qLTsMOHVC04uSz50rDjlx08qO+4tCudjFHrWFBzC0+NWYMHGcL1zoigU41XP/8WJq9Dv+SneFJ7CpKc7fyCLN+3Dte/MxJBXp+HQseBTggje4FpReDBGcdFr00K1AJcMZitG9VBCZfFYeDUE42bW7G9iYuqLAkzWdLSoBBlphIx06dd5gdvxI/3RqqrsgaPh6ZzIm6ec5Ah3qyw41RP6xiPIx3DSI99jxHuzghMgxXE/RqF/V1KvzoqiUEzKDzonGU7rbHFpeQsi6Ko/LT8p5pkmJW5Dl72yKNyMdahGFIWOGWt24aWJqxwfnywTaMKG20HneDitbBUtCnfP4b1p6yT3V0hxU8fGLd6KL+eVnzcVZJX1KoxeFIWOYaNm4KWJqx0fHwaL4khhCZ4asywUqbat4tUYhdPOYtiyz3rB+l2HcO+nCwIdfwkDbh717z6ah28X6nOZqnl3vOo8OUEUhWJUP1wnvZO3f16Dt35eW5ZRMhnwqko4tQSKS8ofp0pvhMl/fd/ni/Dl/M2Yu35P0KIEimqrXVVxYXI9SdSTYlT3jJ28xNHecGFx5e4pAs4q7artBypEnKhq4JnDlwK9sqNcUdjYN1neBVEUIcfJOxxNIxKm3msivKowTm7BQIM058osCjXFKCWJXhNPKFXcn0rF+ymuJ8WofkmcFBcdJ0kmN3vYBrP1qErBEabghCTpzHqO6mdy0/uz8cIPK12XE6JXRRRFPJyGzCn3KzqxKNK8XT/XD1S6elSgKvurXXn8sArD5At3wp/+uxAPfbXY8fGqFcWBY8V45cd8pWUGjSiKODiNelFvUdgvMOrG+Wz2Rnzwyzq1AnmFrnurKtuqqkZAVQNvR57C4tKytZCF+Hw+dxM+mrnB8fFhtbzD1M8TRREHpw2Myt7J5BU7UFhiv7zoGMWuQ4X4y+ilyuTxE1XWkKqnYee55u84iNv+Pdf1OW98fzZOe2KC63LiEe1QXPXWTPz12+R8T1Qga4UkJlBFQUSDiGglEeUT0UiT/S4jIiaiPL9kc9qjVfXKzV2/Bze8PxsLHWSRDMNcDreoGmBUpbftKK4/f7EI45duN/zNjsLRL3/rZnldZsastbvjWjrvTVvnuOywMP3XXdi6z34ivTCNG8USJpdgYIqCiNIBvAZgMIAOAIYTUQeD/WoC+AOAmX7K57RHq+Kdu+Ojebj09V8cHbv/aFGFZVNV9ZgGvTQVH05fp6QsPXrdZrfybt9/1HB7/xem4A8fz3co1XHsjBWYuS3dvB/tH/4e2/YZX2civl6wGVe8OR1fztuMw4XFuOfTBZixZrdzYULI8Ldm4DyDiLVEhHUsL0xiBWlRdAeQz8xrmLkQwCcAhhrs9wSAZwE4qyEOcTyYreDpjlm81fGxpz72Q4WZ5aqyUK7YdgCPfuONi0KfesBO5f16/mac8bdJcX8fXWHmrH3svA5m747b3uv6XYccHhdZwOmD6evw2uR8Jcv1zt+wB7kjx2Du+vAonP0O3vUwNcixhEmuIBVFUwAbY75v0raVQURdATRn5u/MCiKiW4loDhHNKSgoUCJc0K4nlbjtMf1j/Eq8MMF5/isn2FHUM9d6nzDPzvtgtq/VUjbsUrsyXzT8eNGmfXhtsvMV1GKZuiriGvtppZo6FxRGc3i+WbAZnf/6g+fpTczCwsPkEgtywp3RHSq7M0SUBuBFANcnKoiZRwEYBQB5eXlK7q4TRbFg4158v2SbitMrxe0L9+pk70P99C+DnfuvesKU4Tls3ENTRWFB1gnLtuOWD+cYH29ZivJ4MaExOhYWovZMGY9+sxT7jhTh0LFi1KmWFYgMoigibALQPOZ7MwCxPoKaADoB+ElzS5wAYDQRDWFm41qkECe98Item+aBJO5R+b55NYNaX66d++9HhbJzCjPZrQxQLt2yz/rJAoTKJnaGp0FTRfSavMrGak0G+8d4JW6QrqfZANoSUSsiygIwDMDo6I/MvI+ZGzBzLjPnApgBwBclAaiL4w8DKqMnsjP8eWXsWAl+DEbaaQzNxygSH59mUtsTibFt31Gs2n6gwnbV7cctH87BP36IuCPnb9iLmWv8Wy+jpJTx6DdLsHG3WvdclNnrdpfd5yBzMYUpBU9gioKZiwHcCWA8gOUAPmPmpUT0OBENCUquKH64M/xC5ft2tKgUw0ZNV+67dWNRqLq+I4Ul+GzORsPf7HQczKOeEpfjpm3q8fQkw1xVVhq8x79dZvm5Tlh2PPx3+ppduHLUDMsyumXBxr34cPp63P3pAsPf7Xby9OMEl78xHQe1NdKDbKvD1FcNdB4FM49l5nbM3JqZn9K2PcrMow32PccvawJQ30tdtMn+fAhVqH7ZZ6zZjU171C78rq+sdnpTqlwffxu7HPd/vijOOayXY9ZQWSnHrFH/X34BVmzbb10YG7w7bW05BWAXv6zw6NhIcUkp9h4urPB76wfHYp+iGe1B9uqtvNd+ySczs+NQVFKq9CEMeXWawQIn/vDPH1dj6ir7kSlb9x3Bde8Gs1aznVuvqn0qOHAs7m923gXzqCcLFoWJpnht8q8Y9NLPlmWxUmYsbhZsuuLN6b6kjEnXNMXCTfvQ5XHjmesFB+M/Sz1mt8Zr3Wd+biuKQqEwJoii0BF9cANfnIoXFYeE/lpwUGl5Vvlo5gZHDf6rP+bHVTCqXbdGDSgzW5pgpmpCoVmlVed6siOR/7jpHM1dv8eXlDFmYziqCXScwMKp/ZJOFIWO2FfwfcW9ozD5HK2QYZILRHVd1dfHnQeP4cPp69Hj6UlYvtXc1WKlEc8dOQY7DsRXOrkjx2CcSWjzwWPFlv33Zj1BK+2Onw2hnmSIYEq3kKNG1S0Mss5aObdfz0sUhY7YSqra5xqmKAYrpPmYNEp/ay7+1y9lvdNEs5GtVpZV25xbdDd9MAe3/V8k0d9rk/Mx1CQU2nyMworryb58qspMhiAOS4rCYlm/5O8sm7VuhJd1dvLKHabL0DpxPXm1rosoCh3lFIXJg/pu0Rb8okvalogk0xOmFoVqzCpFoh62VX1eWOI8qR4ATFqxAwDw3PiVpska3c7MtnvXl2/dj/1HzQdvrTYgYbYoJq/YgcLiUktXYmVMhplx1dvmKeS8tChueG82Fpi8R5YUhe6N8iqRoCgKHbHvl1mFv/M/8xO+ZHrCXAmNMLMoVPdczO5MIkVhtdd34GgxdtkY5HRKsUm3PJGsew4VJmz09Qx++Wdc+075MahPZ2/AviP2I3/C+orOWLMLN7w/G89PWKmsKbQSuRdkBtdEz2LtzkP4/X/cJ7y0gqyZrSO2UXITAWJEKo1R+LmSX1qC7oxVBXzXJ5G4+3XPXGBVKkeYuW8SidrV4voTuSPH4L7z2uOOvm0AoIKF8+cvFmPyigK8cW039Hv+J/xaYC2ZYFg7M9Gghq17j1pSZlY6D0eLEluYwY5RmJ/8z58vwqx15RMyevX4xKLQEdt5DcNqdUGSbtJCq65AblxPDtZ28hTTFB4KZX1p4ipMM3F/btNSr1tVEkA4OzNHCkvKUudXyUyzpMys7ZP43EEuavTEd8vwzYL4WX7NLFfViEXhIyHtrMXFzKJQ3fP0w/XkB4ncRirvW1EJ42oT96eTYIwwWhQnP/p92eeM9DQcUWQJWLnWlyauxmXdmqFn6/qJC1TMul2HcdcnC1BYXIq+JzVEgxrZ5X5X7fEwQywKH0m2JRfNoktUX4tZnTVTFMUlpaHJy7Vw416c+tgPStKMq8BJQxImpQtUlOc/Mzfgkn8lXtTLyjthRVF8MW8Thr/lLD1JYbGaHv99ny/CHR/Nq7C92MCU9ur5iUWhw8t6EpL2zBILNu7F2p3xXRa+up4MujN3fTIfoxduCZWVtizBfA/A2x77v34qnw6+tJRtNxxhe0eddgLCMKu53cPj8Nxlp+LyvOaJd06AUWCC0b3x6pLEovCRZBqjuOi1afh87qa4vwftevpmgTolocoisRJO7GXj9Oz3K8t9Ly4txaFCeyHBVu6Fn1aH05xrVtz3Tt7h0lLGg18txopt+3G4MPFqepNX7rB9DiOyM9MrbDMao5DBbJ/wsjEPU+/XLdFK9suvOy1VmIT4eHNiG7qPZ21A6wfHui4HADLTE1cnPxvZQ8dK8PpP9hadstJ4+unq89KisFP0vZ8twNqdh7Dz4DH8Z+YGDHrpZ3R4dHzC48Yu3oZez/xYYbvd98Aovb+fz0EUhUs27j5sKcwOCOdAoVOYgS17j+Cqt2bivv8aZ1y1g9k7r/q2Rc+151AhXv3R+ep9sRX17k/mx017HYufb8C2/UeVLXsaZePuw7j0jelKyzTD6YCtFUvETn38ct5m3P/5QkuzwvVs3nukgmIoshmqN2vtbuSOHFP2/deCg4ZlBD7hjojqElFHIjpRW6ZUAND72cm43WCgyYhUUhSlzDigLWS/ekfFhXLsYtbDUt0LL2XGvA170PWJCdi811669Om/Hl+gJ7Yx+nqBtczAYX8HEvVSX5iwynRWumqcBk1YeWfsvlez1+1B979NciSP/r66Wc9l5ppd6Pf8FMN3NxDXExHVJqIHiWgxIivMvQngMwDriei/RNTXG7GCw8mN/smiHzJR2cyMySvMy8odOcZwBbM9hwrx0FeLLcmhgpJSLmv03Caxe2/aWrz189r45/JAUSzd7Gy50dgIGCdh7FNWFpgmJwyaRO2yWSp2L3BsUVgao3BSrhoLx42isDMvRhWJLIPPAWwE0JuZ2zPzWcycx8zNATwDYCgR3eS5lD7i5DXISDR1WCPROzZuyTbc8P7shOUYLTv59Ljl+GjmBktyqKCU1fWO//rtMtOevWpX7JHCEiUGuhMF9vS4FbjcR9eNXRK5LvxWFI4bZkuD8o6KdoT+XIUOFUVJKZsmePTqkkzDY5l5gMlvcwHMVS5R0Di401b9lolM3a0W1l4AjHtZexSt6GWV4aNm4LxOJwDwPi22andNl8cnoFpWxSgSu2zacxgZaWkY8ur/bB1nlq00aBLdaiuT3QBg6qoCTFlVgEcu7OBKHqeKYseBo1iyeR86Na0ddx8/3YD66zCaA2GFRJZI4FFP2hhFdyI6O/rn9uRENIiIVhJRPhGNNPj9XiJaRkSLiGgSEbV0e04nFJeUIj/GD69v8I8UlWDu+t36wyqQ6MW06o81qjwHbCaSc0thSWnZin3b9x/Fb/75P89cKl5ECh22GTZqxKCXfkb/F6YoKSssxGuYmRn/nrHecOlRI657dxbe+V95d+L6XYewbIu9ZVydKoq7PlmAC/9prsB9VRS6czk9t5+zsWOxpCiI6GYAUwGMB/BX7f9jbk5MROkAXgMwGEAHAMOJSN/9mA8gj5lPRcQN9qybczrl2fEr0f+Fqdiw6zC+XbgFrR6oGE556euJ3Qlm70ZxSanll8Co8kQXgw+CXYcKsXjzPoy2OKBrl9ixgPwdB/Cgj2MxlY14DdjSLfvx8NdLsP+o8/esz3M/4fxX7C3j6mUIqJ+uJ30n0Om5i0tKExwbbNTTXQBOB7CemfsC6ArA/iLM5ekOIJ+Z1zBzIYBPAAyN3YGZJzNz1E6fAaCZy3MmxMhHGx0T2HXoGMYt2eq47HjvfGFxKdo8NA5//36FpXKMFMqhY4l7tW9O+dWTaJAoOTWzE+/kgE9mb8CEZdsBAG//vBb/8XEsprIR79HbGXwd+YX7cOkoXvagg3I95e84WC56zg6Jwmq9uiSrKTyOMvNRIgIRZTPzCiJq7/LcTREZKI+yCcAZJvvfBGCc0Q9EdCuAWwGgRYsWLsWqSPQZpxEhy8KkqnjEGyi0O7Bl5KKy8tI/PW4Fdh8qxJGiEjw+tFPc/Zz24mpVzbR9jJVzTVy+AxOX78C6Zy5IKTdPMmE1YAMAPpm9MfFOFvHKoth58Jil9ShUEet66v/CFMflFJWUmg9mB6woNhFRHQBfA5hARHsAuPUzGF2u4WUS0TUA8gD0MfqdmUcBGAUAeXl5rm6V0Y2ONsLpaYQsgxmSbsresf8ollrIERSLUbSN1eHkN6euAQBTReG4F+fgMLuJ08I9C8EeR4tKUFRSippV7CtYv3Ey0QyIWKdWVpuLh1eKIu/JiZ6UG4+xi7ZiQMcT0LROVVflfDV/M35MEELvBVYVxS3MvBfAY0Q0GUBtAN8nOCYRmwDEZstqBgPlQ0T9ATwEoA8z+xubpxF9WQ8dK3YVqmnU6x/88s/YdcjaAKFenljcVEY9TmO8/UhrHfYJa3a44JWf8WvBIc8XUgqSUgbSXbyaYckM7JbHvl2Gf/6Yj7mPxA0ktcRz41ea/h7IzGwi+g0RFQBYRESbiOhMZp7CzKO1cQU3zAbQlohaEVEWgGEARuvO3xWRSX5DmNlTNXqksATPfr/CsDcdbZuuHDXDNFFeIozeebtKAgB2HizElFVuh4gq8vPqAhSVlDoO3Xv7f2vw+k/2UkbYbvhTo90AEMzEqUS8PGl1uVQRUbycz2CG28V5ckeOwRcu6qxKnNR1uwQVHvsUIpPtmgC4FMDTqk7MzMUA7kQkgmo5gM+YeSkRPU5EQ7TdngNQA8B/iWgBEY2OU5xr9hwuxLvTjGcHq+rFbth1CP8Yv9J1uOfHszZgxLuzcCTGX2+306aXYcaaXbj2nVl4eeJqFDmsnDPW7E44ID928VYsj3G12T1VMmXgTWamrCrA6IXHDXynDbbbuqOi7unTrwv2SaQoipl5BQAw80wANVWenJnHMnM7Zm7NzE9p2x5l5tHa5/7M3IiZu2h/Q8xLdE6TOlVx0gm1DH9TpSgWbtqHVyfnKxtEc5PaQt/T26HNuF2765Bji8IKt380D4NfPh4iaffeppDnqYxLX/8FSxymFPGKEe/Owh8+nl/23allMG7JVlezuVW8i1kZxydX5u84YGgxpQpeVY9EiqKhNuntXiK61+B7ShFvvE514xRJJOh+Untx7FiCTZNC72KLRlKlE7l2Fzz53bKES4KWnbcSj1FEmbt+Dx4bvTRoMUxxGuBwz6cLcaOFtDTxUDFGkRUzSDJ9TeKJsclMUK6ntxCxIqJ/+u8pRbwBYS8ap7GLt+Ezl2GEsTHVdl1P+gp4PMGf+5ft7f+txYsTVlna1047cLiwOGUGN/V4nAXFNW7uuz6P19jFW01T85eUMuZv2BP5rKDuxUYqBr3U6wNfqptj4ieJcj391S9BwkC8uqo6e2mU+11OTHIz0KfvIUYbgjQiJYrxmMWwVzsV18pCMcnKgaPFePK7ZUGLYcjSLfvKJjw6IVMX9nT7R/Nw/Zm5eGxIR8P9R01dg79/vwKf3tpDyYS78orCdXGu+HiWujkmRng1hmeqKIioI4DW0TEDInoRkdBYAHiVma0txJAkxOvVuQy88IxY/63d8Ngxi7ZiYMdGaFAjG29M+RX/mhwZ8CthxsTlzhuFKHoFMGHZdsOw2xQ1EBKin728YtsBrNjmfl0PL7jgFXtJD/UYTdYzG6dbu/MgAOBv41YoWfsiduXBVHRdliMg19MzAHbGfD8PwBgAkwE86o1IwUEGNsX3S7baXtzGL9z0th78ajGuf28WAOCZcSvKcvh8s2ALnhyz3LVs45ZsK5dq45YP51RY4Km4pDTQHFVBonL2ctgxmqQaz5L810/5+N/qSJOjaoGk2LXMU71j4tXlJZpw15iZf4n5vp+ZvwAAIvqtRzIFh0Gn/LZ/h9doig5mHysucbQS2JLN+9H9KW9mqO49XIQHv1qMq86In1Lltn/PU2K9COEmwyBKRP+27jtShO8WbcGz35tPKHNCbNUIeowiWUmkKMoNWDNzj5ivDdWLEywhH0+sQHQwu/3DzifJ7/B5IZpYRElUDjLS0yo00PrvD3y5CGMXb/Pk/LHuplTXE14pwkSupy1EVCFRHxH1gPtcT6Ej7JEnetzOWhWERKwpOOi6jKx0qjAmoTeAt+/3rsMSG7GV6mMUQbme/gzgUyJ6H0DUB9MNwAgAV3okU2AYjVGEmfkb9lpeFU8QnHDu884znUZZuGkfej87udw2fYPmZQMeW3aqj1F4RaLw2FmaRXEngOu1zUsB9GDmlPMbJJtF8ZeQT9ICIr05p5lHhdRF7yJxMsZmlajh/fDXi/HvGam9lklgaca1ZHwpF+FkRLIpimTgSFEJho+aEbQYQsjQN2he9vRXbNuPz+ZsTHklAQSUwoOIvtUyyFZImE9EJ2oJ/G70SDbfSTbXUzJwuLAYi0OWx0gIntiJYUeLSjydcb/ncBHu/zw5Z0TbxavB7EQWxS0A7gXwEhHtRmT50yoAcgH8isiku288kSwAxKJQz77D1nI+CZWLqDto2Zb9ttfRFvwn0RjFNgD3A7ifiHIBNAZwBMCqmLWsBSEuu+Pk4JfUz5WbqEWxZItYmyoJKurpuADM6wCs80iOUKBylTghwp7DxorCi4lVQvIQ9TRVyUw331GwR0ApPCoVoibUs0dcT4IRWoMmdS45EEURgxgU6nngy8VBiyCEEAZj18FjeC/OqpKCMwJZMzsWIqpKRO1VnpyIBhHRSiLKJ6KRBr9nE9Gn2u8ztXESzxA9IQj+MHvdHkV0U4IAACAASURBVHR7ciLmbVCT+E+IENTCRQAAIvoNgAUAvte+d3G7fjURpQN4DcBgAB0ADCeiDrrdbgKwh5nbAHgRwN/dnNOCTF4WLwiCkJRYtSgeA9AdwF4AYOYFiITIuqE7gHxmXsPMhQA+ATBUt89QAB9onz8H0I88bM1FTQiCkMwEalEAKGZm1XFsTQHEJuXfpG0z3IeZiwHsA1BfXxAR3UpEc4hoTkFBgWOBxKAQBCGZCXqMYgkRXQUgnYjaEtE/AfyS6KAEGDXL+qu0sg+YeRQz5zFzXk5OjmKRBEEQkoOgLYrfA+gI4BiA/yDSs7/b5bk3AWge870ZKqYuL9uHiDIQWYZ1t8vzxkUsCkEQhIoknHCnDTqPZub+AB5SeO7ZANoSUSsAmwEMA3CVbp/RiKQ0nw7gMgA/sodLVImeEAQhmQkkKSAAMHMJgMNEVFvlibUxhzsBjAewHMBnzLxUSzQ4RNvtHQD1iSgfkZxTFUJoVZJMFsWIni1x01mtghZDEIQQEViacY2jABYT0QQAh6IbmfkPbk7OzGMBjNVtezTm81EAl7s5hx2SKXtsraqZ+OPA9nhv2trQLsbStUUdzJc4+aTn+jNz8f4v64IWQwgQq2MUYwA8AmAqgLkxfylFMlkU0cWAFv5lIC7q0iRgaYzp3KxO0CIILmlZvxru7t/W03NUy5J8T+oIJs145NTMHyTeK/lJKkWhCVuzSiZqV62wXIggKIEAZKR7l+nn1au64scVO/DlvM2enaMyEajriYjWwjgs9UTlEgVIMrme0mKWF00L6VKjRSWlQYsguCSNCJnp3rxfr17VFRee2gQ/rXQ+90nwB6tjFHkxn6sgMm5QT704ARPO9taQjBjlkB5SU6hz8zr4aGbqLz+Z0hCQneGNayhabkj7OUlJYFFPAMDMu2L+NjPzSwDO9UimwEim9zU9VlGEtKaVljJqZkf6ImP+cFbA0giO8DBQokpmpPlJJks+7AS1FCoAgIhOi/mahoiFUdMTiQIkmZICpieB66ladgYm3NsHa3ceQscmSqOrBZ+ZeO/ZuPM/87Fi2wFlZdbQOhFeV7v2jWpi5XZ1coeZoFe4ez7mczEiK91doVyagAlnc2tMrOspI4SK4omhHXHhKY2RlkY4oXaVoMURXNKmYU2cdEJNpYqifvVsAN530MbfczYuf+MXzF63x9PzpDJWo576ei1IGEgig6Kc3zjLw6gUp1zbM9f09990boKLuzbBje/P8UcgwRGxPVTV83Xq1cgC4M8YxWe/7YmZa3dj2KgZ3p8sQIJej+IuIqpFEd4monlENNAbkYIjifQEsjOPP7rMjPApikTk1MhGn3YNgxZDsEGp4laoelZ0MLt8zXt5WBdMuOdsLHxUXRNDRJVife6gXU83MvPLRHQegIYAbgDwHoAfPJIrEJJpjCI7RjlkhtCiSERWRppEu4SUBjWysPNgYYXtKhuhL37Xs6y+3TOgHQ4XluD357ZB9ewM5NSMuKQOHitWeEagJKwpDJIAq4oiWqXPB/AeMy/0cgGhoEimC4p1PSWT3FGyMtKSSjFXJp4Y2gltG9VE/xemlIuiURlR063l8ej6etWz8PwVnSvsozrsu6GmgFIZr6KerHZF5xLRD4goivFEVBNAys2mSqb+RoaLSVB/HdIRLw/rolCa4wzo0Aj/uLxipdeTbdNd9tZ1eZj0xz5OxQo1NatY7a/5Q1oaGVp7pQpq/Me39MC0kdYi61WHfTevV83yuYXyWK2tNyGSufV0Zj4MIBMR95PgkLeuy8PpuXUdH++mt1WnWiYa1vQmEuneAe1wWbdmhr8t/MtAXHVGCwCwPdu3dU51tM6pgVVPDnYtY5h45MIOmHJfX6x75oKgRSkjncjQ2lMxRtE6pzqa1qlqaV8vovmsnlsoj1VF0RPASmbeS0TXAHgYkcWLBIcM6NAIjWo5b6xjK7JdneHlBD2zyl27ambZ73Zn+0ZlzkrCgXszBnZohHrVs4IWoxxpaSiT6eKux5W+CovbznMP6/ygMBN0mvHXAXQmos4A7kdknYgPAaSUL8DDNZHKuHdAO1yt9aqjXH9mLg4XFuOzOZssl+OmDtWq4l0SwURK6M6+bbB5zxFcfJp+eXRz9JExqYIbF6JXpBGhdtVMLH98UNnsaUBN/YiN1hPUE/Sa2cXaynJDAbzMzC8jBWdmG5FTMxsT7jlbWXnN6lZF/RqRQbWqWrje6bn1MPiUxrbKcdPbqlc9y7M5I4kURcNaVfDO9afbVlZhTVPiljBeV1Smqlnp5SxXFf2oMM75SSWCXjP7ABE9AOBaAGO05VFTLre10T3OqZGNto3U6cTYhuHqHi1xebdm6HtSji1/7PDuLdC1+fG1Huw2NfWqZ6FBDW8iQLxq+MLYoKogI+14FZxy3zn419WnmeztD/GstwtOtdeZMSw7BM/x8aEd0apB9UDO3axuVbxgEOEVdqwqiisBHENkPsU2AE0BPOf0pERUj4gmENFq7X+FUV0i6kJE04loKREtIqIrnZ7PDap73rGVsEvzOnju8s6olpVhqyF8+pJTyq0RkG1jIlHbhjXQsGY22jSsUSFRX/N6VR1FhXx8S4+yz7ENn0oqg+upZf3qOKd9ToDSRIh3ry85rRkm3pv83ubreubiwxu7B3LuNCJccppxsIcKArUoNOXwBYBoN3QngK9cnHckgEnM3BbAJBivhX0YwHXM3BHAIAAvEZHvS6apbqDiWQ5uGtjLujXDXf0Sr0LWuVltTLi3T5mS0Sfqq5aZ4SgqpGfr+mWfPdITSi2KMOXG0stSLSv4UFmze92mYQ0fJfGOICyboV2aeG4xBjpGQUS3APgcwJvapqYAvnZx3qEAoqvmfQDgIv0OzLyKmVdrn7cA2AHA0+6WkTZWblHEeUHdNISZ6Wm4Z0C7xDsaXMw9/duhc/M6uP2c1njrujyDg+zhlUVRPVtN+oVXr+qKL28/EwDQqkF1LHv8PCXlOsXoub96VdcAJDmOqja0/8kNcWVecwBAp6a18OKV4XG5BNFXeHlYV3RqmpxZlK3W6jsA9AKwHwC0BtxNop5GzLxVK2trorKIqDuALAC/xvn9ViKaQ0RzCgrUrpZlZfbwoses56SJP//B+4grozPf1b8tvrmjF+4fdBJa1K9mq7wRPVviK63RjeLVWIKqxXPqVcsqU2aZ6eS4B/+3i0/BQ+efXG7btT1aopaFyXP9Tz7+umcaKNYLT22CTk1rOZJLBap62zee1Qp9NFda0zpVy4XaWqV+9Sw0r6d+7kOqroER9GD2MWYuS/5CRBlI0LIR0UQiWmLwN9SOgETUGMD/AbiBmQ3nhjLzKGbOY+a8nBznRofRBUXrzCvDu2JEz5aGx9mJ4InX4fYjDY1q6+iB809G1xblh5fC5NYxIjMjrcw8d+NWbFgzG03rlm/A/nReewzr3iLOEcd5e8TpZZ/DMLirR1XqjHSiMldm3WrO5orMfWQA/n7pqUrkiSWEt10JQScFnEJEDwKoSkQDANwO4FuzA5i5f7zfiGg7ETVm5q2aItgRZ79aAMYAeJiZA8kPHG1MhnRuggbVs/DB9PXlfrfr04+n8X2YwqG8D2XUnoQ9OikzPa0sFYUbWdPTCCgpvy0rPc2XuTheo2pcLj2N0Ll5HbxxzWno3dZ5B86NNRkvus9unrF61bMw75EByB05xrEsvhCwRfFnAAUAFgP4LYCxiMzOdspoACO0zyMAfKPfgYiyEBkw/5CZ/+viXK6IbUuqZJV/Yd+/4XSM/UNvJedRkR7h/kHtTQcbVSfhM+p52m18Ozbx18WSmU4o0e61G0VRp1pmhZn1qTJrXNUwU/T+DurUGNWznQ/S280LFuWZS07B6Dt7Gf5mtyqEqQPw6IUd0L1VvcQ7KiThEyCiNACLmfktZr6cmS/TPru5c88AGEBEqwEM0L6DiPKI6G1tnysAnA3geiJaoP15k8lOw+iSYn2ZNXQve8v61VG7mr3pJF5aFLef0wZPXdTJdTnf/f4sPHphh4T7GfU87botvr6jF969Pv4gerWsdLx0ZfnH/tXtZ+IP57axdZ4oWelpKNFMCjeKon71bHRrWRf/va1n2bb0NLL8HN8ZkVeW9ypsqLIKVZXjVFEM7HgCmsSx+O1aTW6rZx2b7YQZN57VKm76H6+inhKqeWYuJaKFRNSCmTeoOCkz7wLQz2D7HAA3a5//DeDfKs7nhrrVjz9gN72iRKgKOzTzeRdbHAjp1LQ2CksSpwo1Opddn3tmehqqmLgWLjilMS7qWj7dR9cWddGmYQ288mO+rXOVnU+bd9Ksrr3B+1jqa6uznZ5bvmdntZr2O7kR+p3cKO7vQXZgVbmeVJVzYk4NXNK1Kb6cv9nm+eP/ZrdD4+Z5vDysiyvXmxGZcS4u6FxPjQEsJaJZAA5FNzLzEE+kCgj9PX70wg64JCYnUQ1dhIyTahDvOebUzMa6Zy5Ah0e/x+HCkjh7JcZMpiOF1heC8XOkwcwlFu8np41QZkYaOjaojZeHdTFtqM3o2qIOqmUZKzdVFTXI+YWqGnhVeazS0wgvXNnFtqIwe6/s5pxy4xpu27CmksSPt/VpjT/0i1jSfucIs6oo/uqpFCHlxrNalfteTRfLH2th3NanNd6YYhi9awu3j9+sjttRQGar5o3o2RIP6EJD3WAmc7xGy2lbFu2JDe1iLylhLF/dbuz3VkmQ4ZvKXE8Bz6Y3uwzH7qwOjfDDsu22jlFxP9+4phvOaZ9TZg3Hq59eGaKmd4uIqhDR3QAuB3ASgGnMPCX655FMoSX6cGpWycC4u3qXLdkIACMHn2Qp/YX3g2LxX8ojNhSFmZtN9frDZj3YeL3CalkZeHlYl7ghy/EwqrQ9T6xvsKczVPmIVbSxP/3pHPeFuCDo0F9zS9XuaHbk36jr8vDmtd1sHapCUQzqdEK5Onfl6c0N9wtqhbsPAOQhEu00GMDznkgRFizc4/+7qTvG3302Tm5cMVrHyvvguZowkeGQDdeTqpnQVjC7b2a/De3SFI1thicbzRz/+NYeFQbM7XDzWa3QuZnaGbcqItRyA0p8FyXoOTUqTx/rerLbFnsRMn5qszq+LnaVyPXUgZlPAQAiegfALO9FCjdmg1JW3AVeGxRmvfP61a1njK3uY84hszYxUXtppzHq3LwOasSZOX2s2Pm40MMxEWLKxihcHPvBjd1RVOx83VJVvdKgEzmqOn+rBtXxx4HHU+T0btsAPU6shxlrdls63k8XXFAT7orKBGAuVh2HHzbcug2s3J4GNcwHtaxGJsWVIc72fw7viq4trOdUrKrQtZQIs/cqUWW32lsjAr65I/7YQqHFhtWvKuDmPH3auYuwUdXYhH3ypVUm61x41bMz8MmtPS1PvjOalzLroX6YsGw7HvpqiQIJjxNU1FNnItqvfSZEZmbv1z4zMweXkCaEJKoWH918BvJyzSfKHHPREwTiNzC/6dzEVjl++pfNzpRIUVzUpSm+WbAFCzbuNd0vUQWyet/96iVbPcvVZ7TAZd2aoWuLuigpZRTFCWu2E16qqrEJ3vXk7vx39m2jJFrJSGE2rFkFTWqrz2HldJA+EaalMnM6M9fS/moyc0bM55RTEm4riFnPODOd0KtNA3cnsCJDEiY7Mx/MNj+2bvUsfG1iKVjFuqIw/12V28aq9X5FXvOyfFvpafGDDF64sgv6WlzrQj+x1ClBD2a7PX2X5nUqRD46IZ7rSXWfo3bVTHz6256Jd3RAauQcCAmmvnafGnCVL9+8RwbEnS+gkjDct9z61gZ+E8mjyvK32sjZce9Y2ffbO8/CCbWNZ/1GmXpfX3x6aw/TfYDgw2MTKdsnL+qEbi0rrJlWhopUJtWz0uNGEKqwTv26xaIo4tCluf01kkyfWRL4tvXUq55lmIY702Cyzxe/6+l4HQWzCuNXp/T8U07Ad78/K+F+ie5vIoPiEQupUQBvFKSVhukUC9FbLepXwxkWQoq9asTmPzLA0kJdid6da7SliOOh4hnMfKh/XEVh5f6M6NnSdGxt3F29y9b88FJpiKKIIVrJ//KbDvj8NvsmXNBRHoD6Bsaoshll8+zWsh4uPNXeOIil8/ukKYjI0qIybp7xy8O64CarrgyLp7Ejz4k5alen++73Z6GuSQ4jrwZW61bPShgUAlhz37lxeyaiRnaGqRvPSl29d0B7dDbptJ50Qi0M6nSCI/nsIIrCgHrVs8qtSW0VcxeKGhrVMg9x9XKN7yiqs6SaVlalZ3KPG71lp1E32rNZ3aporHML2XGP/HFgO7x3/emJd7RIp6a10aKecb6s7q3qoVZVdYnwAODi2JxfJvfyrn5t8fIwi/NiTEOzzZ/XN3f0wt39jS2bF6/sjB/uOdv0+BIrmtTKK+NDJRFFEYPr8FgFT2zivWfjleHGLpw3rumG7+8yf/nUK4qK21RHVpjPo/BXVSQ6XaLG3uwdsnMpRvte3LUppj9QPpdmInmuPzO3bBA7Mz0NfU+Kv5jkb/ucaF3AMjmNOxKf/ban8vDY5y/vjNVPDU64X8v61SynaDGPuDM/tnPzOjj/lMaGv/Vq3SBu5too0SzGZli5hdFdvJyjFfxK7iHEqXuBTNpPq0W2aVgTTesY99LaNKyBugnC9VTPf4jXEKgkDGMUx89Hpj29CxOEGZtVVjvvldG+Rs8i0f15bEhHS+eb+WC/uKmrzTA6f7GFzMNOSEsjpFnojNlRUGZtq5XnFXcPCyIUlyRu2S3J4ENnSiwKA5zed1WPK547wYpciXoxdjGqdKrWr47iJCmgV5i1MbMe7Icnhpo3vAM7xvcX27kSo8s2ki0rXc2zcHqbjZ6PL8v66r6fnns8eqmOnWVXTWSN51YrJ0ec+2bFu1Bi4UZZeS5+1BCxKGKIPjinE4VMk5DZeJxx464tHGuW9dUJho2Tcosi/m9+xwdEGj7jCly/RnbC3mqfdjlY98wFhrN2bbmeDJ52tFFe87fzsftwIX5csQMt6jtfUyPR+axgpChUra1ihv60/73tTBwtKsHYxVtxdlvr85XiuQoXPjrQ4qJkzudIWMnCYM2iSHwut4hFEUN00lW2Q/eNqslW8RojqybmssfPw0kn1FQii9GLqn72p1nkifNa8PjQjljy1/NsHePlwLqdazGzKNLSCA1qZOOKPOMMok5w6uLTy/nJrT3wiYU5Fm659LRmuEwX2lolMx2XnNbM1n2OV2WtrlwZ36JIjDqLQlxPvlKmKBw2hDWyM3BH39bo1aZijLm9gUznFgUQScGtyrLQh6cSAT0UpuUGzBup1jnOM6CmERnO+TDDzGJw23NzW51VhQobzdB2qpD1irXHifXRoIb15JNOqZKZjn9c3tl1OW7dZPHumpX7acWisKIEUtaiIKJ6RDSBiFZr/+NOjySiWkS0mYhe9VquaGI4p2stEBHuO+8ktG2opjevx1Z4paKXR982rX36gnLrcKgg3uphX91+JobYzFEVSxqR7TEOLyOw3A5mqxqvee+G7vjdOa3LbXNasorZy0HiPhGo806dqqin47JY39cuQT3mkQAmMXNbAJO07/F4AsAUP4SKWhRZiv38gJoBJxVWiV38GEyOl2epa4u6rq4jjezL72XGU7fhsV6trxApW90YRTLhVUipldtiJSjETuSVl+GxQSmKoYgsigTt/0VGOxFRNwCNAPzgh1DRNQmq2FxP1wp+zwdQ1aj4IbfVFN92IbJ/H7xs+OwNZldEpWylereH4zGKJFcULo+P63qycEN/07kJHhh8knn5Vm5vqrqeADRi5q0AoP2vMAuIiNIQWVHvvkSFEdGtRDSHiOYUFBQ4FqqwzKJwF3Jo9HB9tygUnA/wZx6D6iiqWOw2ZN4qCjuuQ+9cT0BF37zz8Fj3sgSKy2543Ptm4b6kpxF+26e1YeLNE2pVwUVdmlh6Z5J6MJuIJhLREoO/oRaLuB3AWGbemGhHZh7FzHnMnJeT43zRluNRT+5ui1cPzm1D4wQ/XAsdm9TGK8O7oqWiUE8j7juvvaX9vGz47BRttK/KRyGupwjuB7Odh8dGMQqH79WmAV4aZi3Jph+PwLN5FMzcP95vRLSdiBoz81Yiagxgh8FuPQH0JqLbAdQAkEVEB5nZbDzDFYUuo56iGPqXFbRAbhsaJ/iVlG9I5yZ4a+oapWVGK7GdtYW9bPjcBiNYmQBmFX1H2vFgdnLribiBFFZxEx4bxaiO2RlkPz5G4d0gRVCup9EARmifRwD4Rr8DM1/NzC2YORfAnwB86KWSAIBHL+yAWlXMMz5awajyqEit4TZqxtk5lRRjiQyboawJcVBceAazy+/8xe96ot/JjZTJom9UnDYxST9G4dlgtvX74vad86MzF9TM7GcAfEZENwHYAOByACCiPAC3MfPNQQh1abdmuNQkP71VjF6SqgoWAHIbNeOEfic1xPwN5suMqiLopTMB43DPU5rWRl5u/AVurGLHJanfs1tL8yV07TL4lMb4YPp6ZKYTiizkHIpH0I9s1oP9cLTIeTCE68HsONdv5764VhTa4SmXFJCZdwHoZ7B9DoAKSoKZ3wfwvueCKcLosTudm5Go3Lj7KqrAt5/TBsO6t0DekxPVFGhCRgiC8o0ssXPa5+CPA62NcZiXbX1frzvqPU6sj3XPXICdB49h9trdjq3o1jk1AGxXK5wNGjpIZBiLW3dN/HkUNiwKgzJq20jRHn1nvUyxJbmevMDgHVGypKitqKfIzpee1gwXdjZOhWyFaLoIP8hUnb7cwTG92zbA+l0blMpRhg2B7jvvJOw4cMxza65BjWwMjpMq2wr3DGiHHifWx5qCg6jn03sSJuI9UjuRfEYWhZ06F1UUbsdbTM/hWcmVGKPehIoxCluuC23Xi7s2Rd/28dcgCBOZiv0YTvznf/lNR/z0p3OUyhHFzrhRm4Y18NXt8ZfADAuZ6Wk4u10Oru/VytUs+qBw27bGcxu5XcvcTvaD6PFeKgqxKDzAqD1Q4nqy0e4dN0d9yPmsCOWD2Q7ITE9DbgPn+aWiDO3SBL8WHMSSzfvLtjm5uqn39XUdri3Ex23jqiL4QV/GNT1aYGgX60o32i6k4szslMbo1bmxV67rcq0sdFImgyaEH2sDqMLJ8rNmBKl2Xh7WFc9fXn45TicWTov61RwtKCRYw239iLckgB30uubJi06xteaLjFEkKXoXw8K/DLQ1OBUPJ70fL2OrVaPa9aQKp7dQr/dCenmVGrcWd7oCK9htEEfUIvGyroui8AB9J8NtA3F5t2Y4rWVdW6vXqe5lPHFRJ7RUOOHLCNUWhSqcNib6DkOSTzlISVyPURg81O659kKZ3c6DSNnw2FRH/9jdTn7LqZmN4d1b2JOh7OVR8/Zc26OlknLMsLt2RCKCbpj1PcVkn5wmVEQ/vvDwBSfj5t4n2izDnQwkUU9Jiq5BcKsonBxeZlEkj+cJvz+3Lc5srW5RpKDbZb1HQdRE+Ihm0T09ty4Gd4q/3nk8lAxmu3xR08sUhWtR4iKKwgMqWBQO7vK/bzoDvW2s/RtPhmRSFE3qVMV/bvF+GU2/EIsi/NTSxg4HdjgBr1/TDTed1QqvXmUtGR+gaDDbtespdVN4pDT6B+fkQZ7VtgHmb9iDn1fvdCRDmevJ0dGCCvQdBBnMDh9Xn9ECaWmEYadH1h9/5MIOto5XkWfJrbLxI6GBKAoPqDiY7exFiDbyztKWe++3DDt+5Ok3Q98ABC2PUJGM9DRfxt/McBv+7IdFIa4nH3Da6Yi28U7eg/NPifhb2zXyZv3usPH+Dafj/kHl8zGpqj9OdW1F15MCYYRQ48S9+PSlp+Cpizs5PqcoiiRFrxjc+qadHH3Jac2w8slBaKVglnEycE77hrj0NPeZf1VSYTBbFEXK4yTKsFaVTFx9hnOrRlxPSYqqQUu3k4HszO4ME+/fcLqj5VHDttqa3qIIm3xCaiCD2UKEStbAnOMwiaGXiw45QSyKykcQkW3iekpSgvaNV1a80hNOn6cMZgt+oCJENxGiKDxAVYNwPOpJsII+VNFN7+7ZS08tS8VQs4qzPF16CydkBo+QIpAPrXggioKI6hHRBCJarf03XGeSiFoQ0Q9EtJyIlhFRrr+SOkNZg6CZFOKysEbFHrxzrji9OU7MiQQC1HKoKPSKSp6j4AWp7HoaCWASM7cFMEn7bsSHAJ5j5pMBdAewwyf5XKH6uYnLwhqqk/DtP1oEAKhVVdVQnjxHQT2p7HoaCuAD7fMHAC7S70BEHQBkMPMEAGDmg8x82D8RnSMNezCoDhPcf6QYgL31i80Qi0LwAj/eq6AURSNm3goA2n+jMJd2APYS0ZdENJ+IniMiw3hPIrqViOYQ0ZyCggIPxbaGqgfXumENAChzgQjmqO5ZPTakI/qf3Ain20wbHQ8JThC8wI9oP8/CY4loIgCjdIwPWSwiA0BvAF0BbADwKYDrAbyj35GZRwEYBQB5eXkpUx2HdG6C3PrVcWqz2kGLkhRUcD25tOzaNKyBt0fkuSpj0WMDccErP2Pj7iOuyhGSgyCMxqSeR8HM/eP9RkTbiagxM28losYwHnvYBGA+M6/RjvkaQA8YKIqwoSqWmojQuXkdJWVVBlQkaFNNrSqZyNR8YuJ6ErzAj9c+KNfTaAAjtM8jAHxjsM9sAHWJKEf7fi6AZT7I5hppD8KBNMyCH0x/4FwM6dxESVkt69tfRdKPSX5Bzcx+BsBnRHQTIm6lywGAiPIA3MbMNzNzCRH9CcAkityJuQDeCkheW7Q/oXIk4gs7oicEP2hcuyrqVc9yXc7ixwYiM6TLAQeiKJh5F4B+BtvnALg55vsEAKf6KJoSerVpgB//2AfnPj8laFEqLed1bOQ4FYhXyGC2YIbTiZ1+ILmePOLEnBpBi1CpefNad4PQgiAcJ5x2jiCkIDJmIiQroigEQRAEU8T1JAiCoIigi1duGgAACMpJREFUrMYqmWm4oVcrz8oXRSEIgpDkrHhisKfli+tJEARBMEUUhSB4TO+2DQCoSy4oCH4jridB8JiHL+yAm3ufiAY1soMWRRAcIRaFIHhMZnoamtezn5pBEMKCWBQe0qlpLeTvOBi0GJWKbi3roltLwwUTBcEzOMWn3Yui8JDvft87aBEqHV/87sygRRAqMak6p1JcT4IgCIpIVbtCFIUgCIJgiigKQRAERYjrSRAEQTAkKyPSlKaHdD0Jt8hgtiAIgkvu6t8OaUS4Iq9Z0KJ4gigKQRAEl9TIzsAD558ctBiekZp2kiAIgqCMQBQFEdUjoglEtFr7bzhDioieJaKlRLSciF4hP1YRFwRBEMoRlEUxEsAkZm4LYJL2vRxEdCaAXoismd0JwOkA+vgppCAIghCcohgK4APt8wcALjLYhwFUAZAFIBtAJoDtvkgnCIIglBGUomjEzFsBQPvfUL8DM08HMBnAVu1vPDMvNyqMiG4lojlENKegoMBDsQVBECofnkU9EdFEACcY/PSQxePbADgZQDTebAIRnc3MU/X7MvMoAKMAIC8vL1Vn0QuCIASCZ4qCmfvH+42IthNRY2beSkSNAeww2O1iADOY+aB2zDgAPQBUUBSCIAiCdwTlehoNYIT2eQSAbwz22QCgDxFlEFEmIgPZhq4nQRAEwTsoiDzqRFQfwGcAWiCiEC5n5t1ElAfgNma+mYjSAfwLwNmIDGx/z8z3Wii7AMB6F+I1ALDTxfHJRGW6VkCuN9WR63VHS2bOMfohEEURZohoDjPnBS2HH1SmawXkelMduV7vkJnZgiAIgimiKARBEARTRFFUZFTQAvhIZbpWQK431ZHr9QgZoxAEQRBMEYtCEARBMEUUhSAIgmCKKAoNIhpERCuJKJ+IKmSzTUaIqDkRTdbStC8loru07YZp3inCK9o9WEREpwV7BfYhonQimk9E32nfWxHRTO1aPyWiLG17tvY9X/s9N0i5nUBEdYjocyJaoT3jnin+bO/R3uMlRPQxEVVJpedLRO8S0Q4iWhKzzfbzJKIR2v6riWiE0bnsIooCkcYFwGsABgPoAGA4EXUIViolFAP4IzOfjEj6kzu064qX5n0wgLba360AXvdfZNfchfIz+P8O4EXtWvcAuEnbfhOAPczcBsCL2n7JxsuITEQ9CUBnRK47JZ8tETUF8AcAeczcCUA6gGFIref7PoBBum22nicR1QPwFwBnAOgO4C/x1vuxBTNX+j8APRHJThv9/gCAB4KWy4Pr/AbAAAArATTWtjUGsFL7/CaA4TH7l+2XDH+IJJCcBOBcAN8BIERmrmbonzOA8QB6ap8ztP0o6Guwca21AKzVy5zCz7YpgI0A6mnP6zsA56Xa8wWQC2CJ0+cJYDiAN2O2l9vP6Z9YFBGiL2GUTdq2lEEzvbsCmIn4ad6T/T68BOB+AKXa9/oA9jJzsfY99nrKrlX7fZ+2f7JwIoACAO9prra3iag6UvTZMvNmAP9AJOXPVkSe11yk7vONYvd5evKcRVFEMFpiNWXihomoBoAvANzNzPvNdjXYlhT3gYguBLCDmefGbjbYlS38lgxkADgNwOvM3BXAIRisFBlDUl+v5j4ZCqAVgCYAqiPiftGTKs83EfGuz5PrFkURYROA5jHfmwHYEpAsStEy734B4CNm/lLbvF1L7w5dmvdkvg+9AAwhonUAPkHE/fQSgDpEFE2nH3s9Zdeq/V4bwG4/BXbJJgCbmHmm9v1zRBRHKj5bAOgPYC0zFzBzEYAvAZyJ1H2+Uew+T0+esyiKCLMBtNUiKLIQGSQbHbBMriEiAvAOgOXM/ELMT/HSvI8GcJ0WUdEDwL6o2Rt2mPkBZm7GzLmIPL8fmflqRFZJvEzbTX+t0XtwmbZ/0vQ4mXkbgI1E1F7b1A/AMqTgs9XYAKAHEVXT3uvo9abk843B7vMcD2AgEdXVrLCB2jZ3BD14E5Y/AOcDWAXgVwAPBS2Poms6CxGzcxGABdrf+Yj4aicBWK39r6ftT4hEf/0KYDEiESaBX4eD6z4HwHfa5xMBzAKQD+C/ALK17VW07/na7ycGLbeD6+wCYI72fL8GUDeVny2AvwJYAWAJgP8DkJ1KzxfAx4iMvxQhYhnc5OR5ArhRu+58ADeokE1SeAiCIAimiOtJEARBMEUUhSAIgmCKKApBEATBFFEUgiAIgimiKARBEARTRFEIQgxEVJ+IFmh/24hos/b5IBH9y6Nz3k1E1zk4LoeIvvdCJkGIRcJjBSEORPQYgIPM/A8Pz5EBYB6A0/h4ziI7x78H4G1mnqZcOEHQEItCECxAROfQ8TUuHiOiD4joByJaR0SXENGzRLSYiL7X0qaAiLoR0RQimktE46OpGHScC2BeVEkQ0U9E9BIR/aKtu9Bd294nxtKZT0Q1teO/BnC193dAqMyIohAEZ7QGcAEiier+DWAyM58C4AiACzRl8U8AlzFzNwDvAnjKoJxeiGRBjaU6M58J4HbtOAD4E4A7mLkLgN7aeYDIzOzeyq5KEAzISLyLIAgGjGPmIiJajMgiOtGxgsWIrCnQHkAnABMiqYmQjkh6Bj2NUX6hJSCSygHMPJWIahFRHQDTALxARB8B+JKZN2n77kAkm6ogeIYoCkFwxjEAYOZSIiri44N9pYjUKwKwlJl7JijnCCJ5iWLRDxwyMz9DRGMQydU1g4j6M/MK7dgjEAQPEdeTIHjDSgA5RNQTiKR7J6KOBvstB9BGt+1K7ZizEMkKuo+IWjPzYmb+OyLuppO0fdshkiRPEDxDLApB8ABmLiSiywC8QkS1EalrLwFYqtt1HCKZUGPZQ0S/ILLc6Y3atruJqC+AEkTSa4/TtvcFMMaDSxCEMiQ8VhAChoi+AnA/M68mop8A/ImZ51g8diqAocy8x0sZhcqNuJ4EIXhGIjKobQsiygHwgigJwWvEohAEQRBMEYtCEARBMEUUhSAIgmCKKApBEATBFFEUgiAIgimiKARBEART/h9eMmVsC9vzJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Convert steps to time in atomman working units\n",
    "steps = output.simulations[0].thermo.Step\n",
    "time = uc.set_in_units(timestep, lammps_unit['time']) * steps\n",
    "\n",
    "# Convert press to atomman working units\n",
    "press = uc.set_in_units(output.simulations[0].thermo.Press, lammps_unit['pressure'])\n",
    "\n",
    "# Plot in ps and GPa\n",
    "plt.plot(uc.get_in_units(steps, 'ps'), uc.get_in_units(press, 'GPa'))\n",
    "plt.xlabel('Time (ps)')\n",
    "plt.ylabel('Pressure (GPa)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8 Reading in dump files <a id='section8'></a>\n",
    "\n",
    "Finally, any dump files generated by the LAMMPS simulation can be loaded into atomman as Systems.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [12.191,  0.000,  0.000]\n",
      "bvect =  [ 0.000, 12.191,  0.000]\n",
      "cvect =  [ 0.000,  0.000, 12.191]\n",
      "origin = [-0.021, -0.021, -0.021]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>atype</th>\n",
       "      <th>pos[0]</th>\n",
       "      <th>pos[1]</th>\n",
       "      <th>pos[2]</th>\n",
       "      <th>atom_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.003687</td>\n",
       "      <td>-0.060490</td>\n",
       "      <td>-0.023425</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.054033</td>\n",
       "      <td>2.032270</td>\n",
       "      <td>1.997684</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2.055068</td>\n",
       "      <td>0.018910</td>\n",
       "      <td>1.980677</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1.997684</td>\n",
       "      <td>2.080987</td>\n",
       "      <td>-0.002075</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>4.081188</td>\n",
       "      <td>0.036592</td>\n",
       "      <td>-0.051009</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>103</td>\n",
       "      <td>1</td>\n",
       "      <td>6.042340</td>\n",
       "      <td>10.090836</td>\n",
       "      <td>8.060340</td>\n",
       "      <td>104</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>104</td>\n",
       "      <td>1</td>\n",
       "      <td>8.037871</td>\n",
       "      <td>8.044527</td>\n",
       "      <td>8.075286</td>\n",
       "      <td>105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>105</td>\n",
       "      <td>1</td>\n",
       "      <td>8.068666</td>\n",
       "      <td>10.122228</td>\n",
       "      <td>10.109220</td>\n",
       "      <td>106</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>106</td>\n",
       "      <td>1</td>\n",
       "      <td>10.135773</td>\n",
       "      <td>8.173499</td>\n",
       "      <td>10.091238</td>\n",
       "      <td>107</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>107</td>\n",
       "      <td>1</td>\n",
       "      <td>10.139345</td>\n",
       "      <td>10.137333</td>\n",
       "      <td>8.031605</td>\n",
       "      <td>108</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>108 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     atype     pos[0]     pos[1]     pos[2]  atom_id\n",
       "0        1   0.003687  -0.060490  -0.023425        1\n",
       "1        1  -0.054033   2.032270   1.997684        2\n",
       "2        1   2.055068   0.018910   1.980677        3\n",
       "3        1   1.997684   2.080987  -0.002075        4\n",
       "4        1   4.081188   0.036592  -0.051009        5\n",
       "..     ...        ...        ...        ...      ...\n",
       "103      1   6.042340  10.090836   8.060340      104\n",
       "104      1   8.037871   8.044527   8.075286      105\n",
       "105      1   8.068666  10.122228  10.109220      106\n",
       "106      1  10.135773   8.173499  10.091238      107\n",
       "107      1  10.139345  10.137333   8.031605      108\n",
       "\n",
       "[108 rows x 5 columns]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = am.load('atom_dump', '100000.dump', symbols='Al')\n",
    "print(system.box)\n",
    "system.atoms_df()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**File Cleanup**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "for fname in Path(pot_id).glob('*'):\n",
    "    fname.unlink()\n",
    "Path(pot_id).rmdir()\n",
    "Path('atom.dat').unlink()\n",
    "Path('nvt.in').unlink()\n",
    "Path('log.lammps').unlink()\n",
    "\n",
    "for fname in Path('.').glob('*.dump'):\n",
    "    fname.unlink()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
